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A.  Research  Objectives 


The  research  objectives  in  this  contract  are  to  develop  and  test  new  integral-equation 
methods  for  solving  electromagnetic  scattering  problems.  One  important  application  of  the 
methods  is  in  calculating  Radar  Cross  Section  (RCS).  Because  the  methods  are  designed 
to  be  extremely  efficient  computationally,  they  will  ultimately  be  able  to  be  applied  to 
calculations  of  the  RCS  of  realistically  large  vehicles  and  platforms.  The  specific  research 
objectives  for  the  program  were: 

(1)  Develop,  test  and  evaluate  codes  to  calculate  electromagnetic  scattering  from  conduct¬ 
ing  bodies  in  two  dimensions,  employing  the  Fast  Multipole  Method  (FMM)  originally 
developed  by  Rokhlin  (Ro85,  Ro90)  for  acoustic  scattering.  The  FMM  reduces  computa¬ 
tion  time  for  electromagnetic  scattering  problems  from  0(n3)  to  0(n4/3),  and  potentially 
0(n  log  n)  per  iteration,  where  n  is  the  number  of  sample  points  on  the  boundary  of 
the  scatterer.  Using  complexification  and  extrapolation,  the  FMM  can  be  made  globally 
0(n4/3),  and  potentially  0(n log  n). 

(2)  Develop  the  theory  of  impedance  boundary  conditions  so  that  ?  mattering  may  be  com¬ 
puted  from  conducting  bodies  coated  with  thin  dielectrics. 


(3)  Develop  methods  to  improve  the  accuracy  of  integral-equation  computations  of  electro¬ 
magnetic  scattering.  Extend  and  improve  a  scattering  code  using  fourth-order-convergent 
quadrature  formulas.  Development  of  this  code  was  begun  under  a  related  prior  IR&D 
project  in  FY1988. 

(4)  Explore  approaches  to  computing  scattering  in  three  dimensions. 
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B.  Research  Accomplishments 

(1)  A  FORTR  \N  code  (FMMl)  implementing  the  Fast  Multipole  Method  for  the  two- 
dimensional  calculation  of  electromagnetic  scattering  from  conductors  was  developed,  tested, 
and  evaluated.  T:\e  code  represents  a  significant  accomplishment  because  of  its  extraordi¬ 
nary  computational  efficiency.  The  operation  count  has  been  reduced  to  0(n4/3)  ( globally ) 
by  using  the  code  in  conjunction  with  the  mechod  of  complexification  and  extrapolation. 
Conventional  iterative  methods  have  operation  counts  of  0(n2)  per  iteration,  and  Gaus¬ 
sian  elimination  is  globally  0(n3).  The  FMM  for  electromagnetic  scattering  is  reported  in 
our  publication  En9l  (see  Section  C  and  Appendix  A). 

(2)  The  theory  of  impedance  boundary  conditions  has  been  developed  for  multiple  thin 
material  layers.  This  theory  allows  for  the  calculation  of  electromagnetic  scattering  from 
coated  conductors. 

(3)  A  code  calculating  electromagnetic  scattering  from  conductors  has  been  developed,  us¬ 
ing  an  integral-equation  formulation  employing  highly  accurate,  fourth-order-convergent 
quadrature  formulas.  This  is  perhaps  the  most  accurate  integral-equation  code  in  exis¬ 
tence  for  electromagnetic  scattering.  Development  of  this  code  (SKIEl)  was  begun  in  a 
related  prior  IR&D  program  in  FY1988  and  reported  in  our  publication  Mu89  (see  section 
C  and  Appendix  A).  However,  the  code  was  extended  and  perfected  under  this  contract, 
because  many  of  its  features  and  subroutines  were  necessary  for  the  proper  design  of  the 
code  implementing  the  Fast  Multipole  Method  (Accomplishment  (1)). 

(4)  The  code  (SKIEl)  as  perfected  and  extended  under  this  contract  was  used  in  an 
IR&D  investigation  of  the  resonance  problem  in  electromagnetic  scattering  computation, 
as  reported  in  our  publication  Mu90  (see  Section  C  and  Appendix  A). 
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(5)  Preliminary  exploratory  work  was  performed  on  the  extremely  difficult  problem  of 
three-dimensional  scattering  computation. 

C.  Publications 

[1]  En91:  Engheta,  N.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990)  “The  Fast 
Multipole  Method  for  Electromagnetic  Scattering  Problems”,  submitted  to  J.  Appl.  Phys.. 

Included  in  Appendix  A 

[2]  En91a:  Engquist,  B.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “Higher 
Order  Impedance  Boundary  Conditions  for  Electromagnetic  Scattering,”  to  be  submitted 
to  Journal  of  Applied  Physics.  Included  in  Appendix  A 

[3]  En91b:  Engquist,  B.  (1991),  “Effective  Boundary  Conditions  for  Helmholtz’s  Equation 
with  Thin  Layers”,  to  be  submitted  to  SIAM  J.  Numer.  Anal. 

[4]  Mu90:  Murphy,  W.  D,,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “Solving  the  Resonance 
Problem  in  Electromagnetic  Scattering  Computation”,  J.  Appl.  Phys.  67  (10),  15  May 
1990,  6061-6065.  Included  in  Appendix  A 

[5]  Mu89:  Murphy,  W.  D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1989),  “Numerical  Second- 
Kind-Integral-Equation  Solutions  of  Electromagnetic  Scattering  Problems”,  Electronics 
Lett.  25,  643.  (This  is  a  relevant  publication  produced  under  a  related  FY1988  IR&D 
project).  Included  in  Appendix  A 

D.  Personnel 

•  WILLIAM  D.  MURPHY,  PhD.  Member  of  the  Technical  Staff,  Computational  Sciences 
Function,  Rockwell  International  Science  Center.  Program  Manager  and  Co-Principal 
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Investigator.  Dr.  Murphy  has  built  on  the  foundations  created  by  Dr.  Rokhlin  for  acoustic 
scattering  to  develop  effective  new  numerical  techniques  for  electromagnetic  scattering. 

•  MARIUS  S.  VASSILIOU,  Ph.D.  Member  of  the  Technical  Staff,  Computational  Sci¬ 
ences  Function,  Rockwell  International  Science  Center.  Co-Principal  Investigator.  Dr. 
Vassiliou’s  role  has  been  primarily  in  the  computer  implementation  of  the  new  numerical 
methods. 


•  BJORN  ENGQUIST,  PhD.  Professor  of  Mathematics,  University  of  California,  Los  An¬ 
geles.  Co-Principal  Investigator.  Professor  Engquist  has  developed  the  theory  of  homog¬ 
enization  and  concentration  allowing  the  calculation  of  scattering  from  conductors  with 
thin  coatings. 

•  VLADIMIR  ROKHLIN,  Ph.D.  Professor,  Dept,  of  Computer  Science,  Yale  University. 
Co-Principal  Investigator.  Dr.  Rokhlin’s  role  has  been  in  mathematical  analysis  and 
program  design. 


E.  Interactions 

(1).  Papers  Presented  at  Conferences: 

[1]  Engheta,  N.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “The  Fast  Multi¬ 
pole  Method  for  Scattering  from  Electrically  Large  Objects,”  accepted  for  presentation  at 
the  1991  National  Radio  Science  Meeting/  IEEE/AP-S  Symposium,  University  of  Western 
Ontario,  London,  Ontario,  Canada,  June  24-28. 

[2]  Engquist,  B.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “New  Methods  for 
Developing  Higher-Order  Impedance  Boundary  Conditions  on  Curved  Surfaces,”  accepted 
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for  presentation  at  the  1991  National  Radio  Science  Meeting/  IEEE/AP-S  Symposium, 
University  of  Western  Ontario,  London,  Ontario,  Canada,  June  24-28. 

[3]  Engheta,  N.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “  What  is  The 
Fast  Multipole  Method  (FMM)  and  How  can  it  Help  You  Solve  Electromagnetic  Scattering 
Problems  More  Effectively?”,  accepted  for  presentation  at  the  Progress  in  Electromagnetics 
Research  Symposium  (PIERS),  Cambridge,  Massachusetts,  July  1-5. 

[4]  Murphy,  W.D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “The  Fast  Multipole  Method 
(FMM)  for  Electromagnetic  Scattering  Problems”,  presented  at  the  URSI/National  Radio 
Science  Meeting,  Jan.  3-5,  Boulder,  Colorado. 

[5]  Murphy,  W.D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “Solving  the  Resonance  Prob¬ 
lem  in  electromagnetic  Scattering  Computation”  ,  presented  at  the  URSI/National  Radio 
Science  Meeting,  Jan.  3-5,  Boulder,  Colorado. 

(2).  Other  Laboratories  and  Agencies 

Harold  Brooks  of  the  Naval  Weapons  Laboratory,  China  Lake,  CA,  has  been  made  aware 
of  our  work  and  is  interested  in  receiving  copies  of  the  code  FMMl. 
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F.  New  Discoveries 

W<  have  developed  the  first  Fast  Multipole  Code  that  solves  electromagnetic  scattering 
problems.  Previous  implementations  of  the  method  solved  acoustic  or  potential  problems. 
We  have  developed  a  new  method  of  calculating  electromagnetic  scattering  and  RCS  that  is 
both  highly  efficient  and  highly  accurate.  The  Fast  Multiple  Method,  used  in  conjunction 
with  complexification  and  extrapolation,  reduces  computation  time  for  electromagnetic 
scattering  problems  from  0(n3)  to  0(n4/3),  and  potentially  O(nlogn),  where  n  is  the 
number  of  sample  points  on  the  boundary  of  the  scatterer. 

We  have  also  developed  new  theory  on  impedance  boundary  conditions  for  multilayered 
dielectrics,  and  begun  work  on  the  necessary  theory  to  better  handle  corner  singularities. 
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G.  Background  Information,  Theory,  and  Computational  Results 

Gl.  The  Fast  Multipole  Method  (FMM)  Applied  to  Electromagnetic  Scattering 

Our  paper  En91  (reproduced  in  Appendix  A)  should  be  co  '.suited  for  full  details  and 
results.  Below,  we  present  some  background  information,  partially  excerpted  from  that 
paper.  Essentially,  the  Fast  Multipole  Method  that  we  have  developed  and  applied  is  one  of 
the  most  efficient  methods  available  for  computing  electromagnetic  scattering.  The  method 
reduces  the  operation  count  for  solving  the  Magnetic-Field  Integral  Equation  (MFIE)  from 
0(n3)  for  Gaussian  elimination  to  0(n4/3)  per  conjugate-gradient  iteration.  If  our  method 
of  “complexification”  and  extrapolation  of  the  wavenumber  (mentioned  below  and  de¬ 
scribed  in  detail  in  the  paper  in  Appendix  A)  is  applied,  then  the  condition  number  of 
the  discrete  system  is  bounded;  consequently,  the  operation  count  of  the  entire  FMM  (all 
iterations,  )  becomes  globally  0(n4/3). 

The  Fast  Multipole  Method  (FMM)  was  developed  by  Rokhlin  (Ro85,  Ro90)  to  solve  acous¬ 
tic  scattering  problems  very  efficiently.  We  have  modified  and  adapted  it  to  the  second- 
kind-integral-equation  formulation  of  electromagnetic  scattering  problems  in  two  dimen¬ 
sions.  The  present  implementation  treats  the  Dirichlet  (TM)  problem  for  two-dimensional 
closed  conducting  objects  of  arbitrary  geometry. 

Consider  a  scatterer  with  n  nodes  on  its  boundary.  Divide  the  boundary  into  p  segments, 
each  containing  n/p  nodes.  Instead  of  calculating  n2  interactions  among  n  current  sources 
on  the  boundary,  consider  each  segment  to  be  a  cluster  of  n/p  sources.  For  segments  that 
are  close  together,  the  exact  interactions  must  be  calculated.  For  segments  sufficiently 
far  apart,  however,  we  may  combine  the  sources  in  each  segment,  approximating  their 
radiation  fields  by  the  first  N  multipoles.  We  describe  each  segment  via  an  equivalent 


7 


Rockwell  International 

Sctonc*  Center 


SC71004.FR 

multipole  located  at  the  segment’s  center.  We  can  calculate  the  contribution  of  each  such 
equivalent  expansion  to  the  field  at  the  center  of  any  sufficiently  separated  receiver  segment, 
and  in  that  receiver  segment  we  can  use  a  partial  wave  expansion  to  obtain  the  field  at 
all  the  individual  n/p  nodes.  The  radiation  field  at  any  particular  node  on  the  scatterer 
boundary  is  the  sum  of  contributions  of  N  multipole  expansions  of  each  of  the  far-away 
segments,  and  the  direct  contribution  of  very  close  segments.  The  method  reduces  the 
operation  count  for  solving  the  Magnetic- Field  Integral  Equation  (MFIE)  from  0(n3)  for 
Gaussian  elimination  to  0(n4/3)  per  conjugate-gradient  iteration. 

We  also  present  a  simple  technique  for  accelerating  convergence  of  the  iterative  method: 
“complexifying”  k,  the  wave  number.  This  has  the  effect  of  bounding  the  condition  number 
of  the  discrete  system;  consequently,  the  operation  count  of  the  entire  FMM  (all  iterations) 
becomes  0(n4/3).  We  present  computational  results  for  moderate  values  of  fca,  where  a 
is  the  characteristic  size  of  the  scatterer.  Methods  and  results  are  described  in  detail  in 
En91,  included  in  Appendix  A. 

Consider  a  two-dimensional  conducting  body  whose  axis  is  aligned  with  the  z  coordinate 
axis.  A  monochromatic  electromagnetic  wave  incident  on  this  structure  with  an  electric 
field  vector  parallel  to  the  axis  of  the  body  is  referred  to  as  the  transverse-magnetic  (TM) 
case.  The  incident  and  scattered  fields  both  satisfy  the  following  Helmholtz  equation 

A  Ez  +  k2Ex  =  0.  (1) 

The  boundary  condition  for  equation  (1)  is  that  the  total  E  field  vanish  on  the  surface  of 
this  conductor,  i.e., 

£‘ot  =  0  on  C  (2) 
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or  more  explicitly 


Einc  +  E»ca  t  =  o  on  C. 


(3) 


The  above  Helmholtz  equation  can  be  reduced  to  the  following  second-kind  integral  equa¬ 
tion 

tf(r)  +  2  /  ^  =  —  2Ei„nc(r)  (4a) 

J  C 

Er\r)  =  J f  d1'  (4&) 

where  r  and  r'  are  both  on  the  boundary  C,  and  G  is  the  free-space  Green’s  function  in 
two  dimensions,  i.e., 

G(k\r  -  r'|)  =  iH<‘>(k|r  -  r')/4  (5) 

with  —  r'|)  being  the  Hankel  fuction  of  the  first  kind  of  order  zero.  See  Co83  for  a 

derivation  of  equation  (4).  Equation  (4),  the  TM  case,  is  often  referred  to  as  the  Dirichlet 
problem  in  the  mathematical  literature  (Kr89,  Co83).  If  we  discretize  the  boundary  into 
n  points,  then  the  above  integral  equation  (4a)  is  converted  to  the  following  linear  system 
(via  Nystrom’s  method  (Kr89)): 

n 

*(r,)  +  2^Aii*(ri)  =  -2I?,"'(ri)  (6) 

j=l 

where  the  matrix  A  =  ( Aij )  isnxn  and  the  vectors  (^(r;))  and  (Eznc(rj))  are  col¬ 
umn  vectors  having  n  rows.  Applying  normal  matrix  multiplication,  -4$  requires  0(n2) 
operations.  The  FMM  algorithm  reduces  this  to  0(n4/3)  or  ultimately  to  0(n  log  n). 
Consider  n  nodes  on  the  boundary  of  the  scatterer.  Divide  the  boundary  into  p  equal 
segments,  where  2  <  p  <  n.  In  each  segment,  there  are  n/p  nodes.  If  the  length  of  the 
boundary  is  £,  each  segment  has  length  L/p.  The  center  of  each  segment  is  located  at 
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Zi(i  —  1,2, ...  ,p).  In  scattering  problems,  each  node  can  be  treated  conceptually  as  if  it 
were  a  source  of  radiation. 

If  we  have  sources  within  a  finite  region  of  space,  the  radiation  emitted  from  these  sources 
in  the  far  zone  can  be  approximated  using  a  collection  of  multipoles  located  at  the  center 
of  the  region  (Ro85,  Ro90,  Gr87).  The  multipole  approximation  converges  rapidly  outside 
any  circle  D  containing  all  sources  and  separated  from  D  by  at  least  one  wavelength. 
In  fact,  once  a  sufficiently  large  number  of  multipoles  is  included,  the  accuracy  of  the 
approximation  increases  superalgebraically  (faster  than  any  negative  power  of  N  (Ro90)). 

Consider  each  segment  on  the  boundary  as  a  cluster  of  n/p  sources.  The  sources  in  each 
segment  axe  treated  as  a  single  aggregate  source,  and  the  radiation  field  of  that  equivalent 
source  is  approximated  using  the  first  N  multipoles  located  at  the  center  of  the  segment. 
For  each  pair  of  sufficiently-separated  segments,  the  radiation  of  the  N  multipoles  of  one 
segment  can  be  represented  as  an  analytical  partial  field  expansion  around  the  center  of  the 
other  segment.  Then  from  this  information,  the  field  at  the  other  nodes  of  that  segment 
can  be  evaluated  using  equation  (16).  For  nearby  segments,  the  direct  contribution  must 
be  calculated  to  evaluate  the  radiation  field.  The  radiation  field  at  any  particular  node 
on  the  boundary  is  the  sum  of  the  contributions  of  N  multipoles  of  each  of  the  far-away 
segments  and  the  direct  contribution  of  the  nearby  segments.  Ro90  considers  the  precise 
mathematical  description  of  the  process. 

To  illustrate  the  above  verbal  description  mathematically,  let  us  consider  the  field  of  a 
two-dimensional  magnetic  surface  current  distributed  over  a  two-dimensional  body.  Since 
we  have  an  exterior  Dirichlet  problem,  a  double-layer  potential  is  needed.  For  the  two- 
dimensional  TM  problem,  a  magnetic  current  density  at  any  point  of  the  boundary  is  a 
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vector  in  the  x-y  plane  tangent  to  the  body.  The  electric  field  of  such  a  current  distribution 
is  given  by 

E““(fl  =  VxjfG(k|?-?'|)K(?')dl'  ‘  (7) 

where  K (p1)  is  the  magnetic  surface  current  density,  p  =  (/>,  0),  and  C  is  the  contour  of 
integration.  In  two  dimension?  equation(7)  can  be  written  as 

Er'(p)  =  (8) 

which  is  identical  with  equation(4b)  if  |K(p')|  is  taken  to  be  'f(p').  Substituting  equa¬ 
tions)  into  (8),  we  get 

gr 1 \?) = m  Jc  — yi)  ik(p')|  ar  (9) 

The  Hankel  fuction  —  p'\)  can  be  expanded  in  terms  of  higher  order  Hankel 

functions: 

Bp'w-F l)=  E  (io) 

m=— oo 

Substituting  equation  (10)  into  equation  (9)  yields 

gr'(ft«)=  E  -P)  <*V(i  m0)  |K(P',  fl-)|  di-  (11) 


This  can  be  regarded  as  the  multipole  expansion  of  the  source  K(p',6').  For  a  discretized 
source  at  n  points  located  at  x'j  —  (p'j,  O'j)  (j  =  1,2,...,  n)  over  a  segment  of  the  boundary, 
equation(ll)  reduces  to 


Er'(p,t>)=  E  ^(i/4)££=M^L^ff«)(^)exp(<m«)|K(X|)|Ali'  (12) 


m=— oo  jf=l 
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where  A /'■  is  the  discretized  element  of  axe  length  containing  the  source  |K(x-)|.  For  a 
given  accuracy,  we  can  truncate  the  infinite  sum  in  equation(12)  at  iV,  and  thus  calculate 
the  first  N  multipoles  of  the  source. 

Historically,  the  FMM  algorithm  was  first  applied  to  Poisson’s  equation  for  n  point  charge 
lines  at  locations  aq(i  =  1,2, ...,n)  with  strengths  /q.  This  is  mathematically  equivalent 
to  solving  the  equation 

n 

A<j>  =  ^  6(x  —  Xi)Ki  (13) 

i=l 

where  6(x)  is  the  Dirac  delta  function  and  x  and  X{  are  two-dimensional  points.  The 
solution  to  equation  (13)  is 

n 

<K x )  -  Ki  los(ix  ~  *«■!)/ (27r)  (u) 

t=i 

If  we  evaluate  equation  (14)  at  each  point  Xi(i  =  1,2, ...,n),  then  this  computation  re¬ 
quires  0(n2)  operations.  However,  if  large  numbers  of  particles  are  combined  into  single 
computational  elements,  then  this  operation  count  can  be  reduced  if  an  approximate  an¬ 
swer  (to  a  specified  accuracy)  is  desired.  When  a  cluster  of  particles  is  “far  away”  from  a 
particular  point,  then  the  potential  of  the  cluster  is  approximated  by  the  potential  induced 
by  a  single  computational  element  located  inside  the  cluster  (Gr87).  In  the  FMM  algo¬ 
rithm  the  computational  element  is  a  Laurent  expansion  centered  at  a  circle  containing 
the  cluster  of  particles.  Given  a  cluster  of  charges  located  at  points  Zi(i  =  1,2, ...,nc), 
the  expansion  is  given  by 

=  Re(J2  log(z  -  Zi ))  «  Re(a0  log(z  -  z0)  +  ak/(z  -  z0)k )  (15) 

«=i  *=i 

Here  p  is  the  order  of  the  multipole  and  the  aks  are  coefficients  chosen  so  that  the  truncated 
series  is  an  accurate  approximation  of  the  potential.  For  equation  (15),  the  computational 
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effort  is  only  O(p)  operations  and  is  much  less  than  0(nc)  for  the  direct  approach.  The 
region  must  now  be  organized  into  well-separated  points  and  very-near  points.  For  near 
points,  the  direct  evaluation  of  (14)  is  used.  See  Ro85,  Ro90,  ar.d  Gr87  for  details  of  decom¬ 
posing  the  regions  into  boxes  of  different  sizes.  When  applying  the  FMM  to  Helmholtz’s 
equation  instead  of  Poisson’s  equation  the  expansion  (15)  is  replaced  by  the  standard  Han- 
kel  function  expansion,  which  is  then  truncated  to  obtain  a  given  accuracy  requirement. 
That  is, 

Eacat(0  0)  =  {  amHm\kp )  exp (*m0),  if  p  >  a  (16a) 

'  '  \'L'!n=-ocPmJm(kp)exp(im9),  if  p  <  a  (16b) 

depending  on  whether  the  calculations  are  to  be  done  outside  or  inside  the  circle  of  radius 
a  (Ro90).  Obviously,  any  prescribed  accuracy  in  the  series  can  be  guaranteed  by  taking 
more  terms  in  the  expansion  at  the  expense  of  more  CPU  time.  We  refer  the  reader  to 
En91  (included  in  Appendix  A)  for  more  details. 
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G2.  Impedance  Boundary  Conditions  for  Thin  Dielectrics 

Higher  order  impedance  boundary  conditions  for  thin  coatings  about  closed  conductors  in 
two  dimensions  are  derived  using  Fourier  integral  techniques.  Using  a  single-layer  potential 
and  these  impedance  boundary  conditions,  second-kind  weakly  singular  integral  equations 
are  derived  that  model  TE  electromagnetic  scattering  problems.  These  integral  equations 
are  solved  using  Nystrom’s  method  and  approximately  fourth-order  convergent  quadrature 
formulas. 

Consider  a  two-dimensional  closed  perfect  electrical  conductor  coated  with  a  thin  layer 
of  dielectric  and/or  magnetic  material.  The  classical  way  of  solving  the  electromagnetic 
scattering  problem  from  such  an  object  is  to  develop  an  integral  equation  in  which  the 
contour  of  integration  contains  both  the  conductor  and  the  outer  surface  of  the  dielectric. 
The  difficulty  with  this  approach  is  that  as  the  thickness  of  the  dielectric  layer  approaches 
zero,  an  ill-conditioned  equation  may  result.  In  addition,  the  size  of  the  discrete  linear 
system  is  twice  as  large  as  in  the  method  described  below.  Our  procedure  will  translate 
the  boundary  condition  on  the  surface  of  the  conductor  to  the  dielectric-air  interface  by 
developing  an  impedance  boundary  condition  on  the  interface.  The  resulting  integral 
equation  will  only  have  to  be  integrated  along  the  interface,  thereby  reducing  the  number 
of  unknowns  for  the  discrete  problem  by  a  factor  of  two  and  possibly  removing  the  ill- 
conditioning  caused  by  grid  points  on  the  conductor  and  dielectric  being  too  close  together. 
Other  work  in  this  area  can  be  found  in  Ro89,  Se89,  Vo89,  Ha75,  Ka65,  Se81,  Se62,  Se91, 
and  Ba90. 
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G.2.1  Derivation  of  Higher  Order  Impedance  Boundary  Conditions 


The  following  treatment  is  taken  from  our  paper  En91a,  in  preparation  for  submission  to 
the  Journal  Applied  Physics.  A  draft  of  the  paper  is  included  in  Appendix  A. 

Consider  Helmholtz’s  equation  written  in  co-ordinates  normal  (n)  and  tangential  ( t )  to  an 
infinite  flat  scatterer  along  the  t  axis,  i.e., 


Q2utot  d2  um 

dn 2  +  dt 2 


+  k\u*°* 


=  0 


(1) 


where  k\  is  the  wave  number  in  the  dielectric  (ki  =  kairy/eijli).  Here  the  superscript  tot 
denotes  the  total  field.  Thus 

u*°*=u  +  uinc  (2) 


where  inc  denotes  the  incident  field.  No  superscript  is  used  for  the  scattered  field.  The 
boundary  condition  on  the  conductor  for  the  TE  polarization  case  is 


dutot 

dn 


=  0 


(3) 


and  on  the  dielectric-air  interface  we  have 


6x1*2*  _  e2  du\° * 
dn  ei  dn 


where  cj  is  the  dielectric  constant  in  the  thin  layer  and  e2  is  that  in  air.  We  assume  the 
layer  thickness  is  6. 

Taking  Fourier  integral  transforms  along  the  scatterer  and  assuming  functions  decay  at 
ioo,  we  can  derive  the  following  equation,  where  '  denotes  Fourier  transformation: 


d2utot 
dn 2 


(  'l  7  2\  *  tot  A  tot 

=  —  k\)u  =  au 


0  <  n  <  6 


(5) 
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Let  Uq0<  denote  the  unknown  value  of  utot  at  n  =  0  (on  the  metal-dielectric  interface). 
Then  we  obtain 


u‘°‘(0)  =  ui,ot 

(6a) 

^“(O)  -  arfl« 

dn*'  -  0 

(66) 

3i‘°'(0)  d™u'°<(0)  „  .  ,  „ 

an  =  9n»+'  =0’  r  =  1'2>--- 

(6c) 

where  we  have  used  equations  (3)  and  (5).  Expanding  in  a  Taylor  series 

in  6,  we  have 

utot(6  —  0)  =  Uq0*  +  ■  H - 

(7) 

afi,7'0)=^'+^'+- 

on  6 

(8) 

Substituting  (8)  into  (4)  gives  (at  n  =  S  -1-  0) 

du^/dn  =  e2%tiJot  +  ( aS)2ulot/6  +  •  •  •)/«! 

=  e2Sa  ut20t/e1  +  0(62) 

(9) 

where  we  have  used  equation  (7)  and  the  continuity  of  utot.  Note  that  the  subscript  2 

represents  the  point  n  =  8 +0.  Taking  inverse  transforms  of  equation  (9)  yields  the  desired 

impedance  boundary  condition 

du\ot  _  -6e2{d2u\otldt2 +  k\u\ot)  [  2 

dn  '  ei  1  ' 

(10) 

In  terms  of  scattered  and  incident  fields  equation  (10)  may  be  re-written  as 

gn  +  %/ei)  [^jjr  +  *?«»]  = 


a“2  ■  -6{e2/tl)[^l+kiur} 


dn 


(11) 
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where  we  have  dropped  the  0(6 2)  term  in  equation  (11). 

Although  equation(ll)  has  been  derived  for  a  flat  surface,  it  remains  true  to  0(62)  for  a 
curved  surface  as  well.  This  can  be  seen  by  making  a  rotation  of  angle  a  between  the  flat 
boundary  and  the  curved  boundary,  i.e., 


t  —  y  cos  a  —  x  sin  a 


n  =  x  cos  a  +  y  sin  a. 


Then  to  leading  order  a  =  y/R  where  R  is  the  radius  of  curvature.  Specializing  to  the 
point  (x,j/)  =  (0,0),  the  origin,  we  have 


or 


and 


or 


d2t 

Oy 2 


=  0 


(z,y)=(0,0) 


d2n 


2 

R' 


l(x,y)=(0,0) 

Using  the  above  and  the  chain  rule  in  equation(lO)  leads  to 

du20t  __  Cf  ,  s'd2U20t  ,  2  5u|ot  2  M.  rf2. 

-sr--^‘>Mi-§r+jsr+kiut 


or 
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Finally, 


for  a  curved  surface. 


To  develop  a  third  order  impedance  boundary  condition  on  a  curved  surface,  start  with 
Helmholtz’s  equation  in  normal  and  tangential  coordinates  as 

3n*  +2H  dn  +  flt!  +*■“  0  (12) 

where  //  is  the  curvature  at  the  boundary  point.  Let  x  and  y  be  points  separated  by  a 
distance  6  on  the  dielectric  and  conductor,  respectively.  Then  expanding  in  Taylor  series 
along  the  normal  yields 

«'*(»)  =  «*“'(*)  -  +  o(s3).  (13) 

Integrate  equation  (12)  along  the  normal  and  apply  equation  (3): 


dutot  ozr  f*  .  ,  &  f* 
dn  +2Hi  dn  ds+dtiJn 


da  +  J  utot  da  +  k\  J  utot  ds  =  0 


where  H  is  treated  as  a  constant  to  0(S3). 
JYom  equation  (13),  we  have 


r 6  di 

Jo  £ 


Inserting  these  expansions  into  (14)  gives 


M0t(*)  ,  ~II(Sdut0t(*)  62  d2utot(x)  d2Utot(x)  2  tot  3 

^T~+2H(6~^^~^'^l^)  +  6~dF~+6klU  (*)  +  °(6)  =  °.  (17) 
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Note  that  since  dutot/dn  =  0(6),  we  have  that 

/  utot  ds  =  6utot(x)  +  0(6Z) 

Jo 

and  from  equation  (12) 

Using  equation  (19)  in  equation  (17)  leads  to 


(18) 


(19) 


+  2H6—^X).  +  (s  +  +  k\(6  +  62H)utot(x)  +  0(6 3)  =  0.  (20) 


dn 


dt2 


The  last  expression  may  be  re-written  as 

Olfll  -  6  (d2u\ot  2  tot  3 

&T“  TT6h{  dt2  +klUl  )  +  °{b  y 


(21) 


The  interface  condition  (4)  may  be  used  as  before  to  obtain  the  appropriate  condition  at 
the  dielectric-air  interface. 

Now  introduce  the  single-layer  potential 


«(i)  =  J^(x,y)(Ky)ds(y) 


(22) 


where  x  and  y  are  two-dimensional  points,  C  represents  the  outer  contour  (around  the 
dielectric),  and  $  is  the  two-dimensional  free  stream  Green’s  function  given  by 


«(*,!,)  =  iff  (23) 

Here  ffg  denotes  the  first  kind  Hankel  function  of  order  zero,  and  k  =  kajr.  If  equation 
(22)  is  substituted  into  (11)  and  the  appropriate  jump  conditions  (Co83)  are  enforced  at 
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the  dielectric  surface,  we  can  derive  the  following  weakly  singular  second-kind  integral 
equation: 

fl Kx )  -2  /  [3$(z,y)/dn(a:)  +  6(e2/el)(d2^(x,y)/dt2(x)  +  fc?$(x,y))]^(y)ds(y) 

J  c 

=  -2g(x)  (24) 

See  Co83  for  more  details  on  derivations  of  integral  equations  in  the  form  of  equation 
(24).  If  6  is  set  to  zero,  equation  (24)  reduces  to  the  TE  polarization  case  for  scattering 
from  a  perfect  electrical  conductor.  If  the  term  d2$/dt2  is  removed  from  equation  (24), 
we  have  the  standard  first  order  impedance  boundary  condition  integral  equation  (Co83). 
Finally,  as  written,  we  call  equation  (24)  the  second  order  (0(62))  impedance  boundary 
condition  integral  equation  for  modeling  a  thin  coating  about  a  metal  conductor.  The 
second  derivative  term  obviously  models  the  curvature  of  the  scatterer. 

The  above  derivation  can  easily  be  extended  to  multiple  dielectric  surfaces.  Suppose  we 
have  m  overlapping  coatings  with  physical  parameters  (6j,e,-, Jb;)  ( i  =  1,2, ...,m).  Let 

6  =  6i  +  H - f-  6m  and  set  =  Cj+i/ej  where  em+i  =  eoir.  Then  following  the  same 

steps  as  above,  we  can  derive  the  following  impedance  boundary  condition  at  the  m-th 
dielectric-air  interface: 

m 

du/dn  =  (d2u/dt2  +  k)u)+g  +  0{62)  (25) 

j= 1 

where  we  have  dropped  the  subscript  2,  and  g  is  now  defined  to  be  equal  to  the  sum  in 
equation  (25)  for  the  incident  field.  An  integral  equation  analogous  to  equation  (24)  can 
be  easily  written. 

In  addition,  fourth  and  higher  order  impedance  boundary  conditions  may  be  developed  by 
allowing  more  terms  in  the  Taylor  series  expansions  (7)  and  (8).  For  example,  the  fourth 
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order  impedance  boundary  condition  is  given  by 

du/dn  +  (c2M)(*3/6)[d2/&2  +  k2][d2u/dt2  +  k2u]  =  g  (26) 

where  —g  now  takes  the  form  of  the  second  term  in  equation  (26)  for  the  incident  field  and 
the  subscript  2  has  been  dropped. 

G2.2  Examples 

Assume  the  incident  field  is  in  the  form  of  a  plane  wave  with  the  incident  angle  /3,  i.e., 

u,nc(x)  =  exp[ifc(xj  cos  +  x2  sin  /?)]  (27) 


where  x  =  (xj,x2). 

Fig.  1  shows  the  results  (expressed  as  Radar  Cross  Section  per  incident  wavelength,  in 
dB)  for  first  and  second-order  impedance  boundary  conditions,  for  a  perfectly  conducting 
elliptical  cylinder  (axes  a  =  2,  b  =  1;  kb  =  5)  having  a  thin  coating  of  thickness  6  =  0.005. 
The  coating  has  material  parameters  e  =  fi  =  14  +  1.4i.  The  figure  also  shows  a  perfectly 
conducting  elliptical  cylinder  without  a  coating,  for  comparison. 
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G3.  Achieving  High  Accuracy  in  Second-Kind  Integral  Equation  Solutions 

In  this  section  we  describe  the  application  of  highly  accurate  quadrature  formulas  to  the 
computation  of  electromagnetic  scattering  via  second-kind  integral  equations.  Our  code 
using  this  method  is  perhaps  the  most  accurate  integral-equation  method  in  existence  for 
computing  electromagnetic  scattering.  The  method  handles  the  singularity  in  the  kernel 
(Green’s  function)  by  employing  quadrature  formulas  with  endpoint  corrections  (Ro85a, 
Mu89)  that  are  0(l/n* )  where  n  is  the  number  of  grid  points  and  fc  is  a  positive  integer. 
Thus,  these  quadrature  formulas  are  convergent  even  for  singularities  of  the  form  logx  and 
xa(— 1  <  a  <  0).  In  our  computer  program  SKIE1,  k  is  approximately  4  (but  can  be  made 
higher). 

In  two  dimensions,  electromagnetic  scattering  from  a  closed  conducting  object  defined  by 
the  curve  7  can  be  described  by  solving  Helmholtz’s  equation 

V2<t>  +  k2<f>  =  0  (1) 

where  k  is  the  wave  number  of  the  incident  radiation.  In  the  case  of  TM  polarization,  <f> 
represents  Ex ,  the  2-component  of  the  electric  field,  and  the  boundary  condition  on  7  is 
given  by 

+  =  f=  -ET  on  7  (2) 

where  E'znc  is  the  2-component  of  the  incident  electric  field.  We  refer  to  problem  (1-2)  as 
an  exterior  Dirichlet  problem.  For  TE  polarization,  $  represents  Hz,  the  2-component  of 
the  magnetic  field,  and  the  boundary  condition  on  7  is  given  by 

0H[nc 

dv 


d<f> 

du~9~ 
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where  u  is  in  the  outward  normal  direction  from  7  and  H'*c  is  the  ^-component  of  the 
incident  magnetic  field.  We  refer  to  the  problem  defined  by  equations  (1)  and  (3)  as  the 
exterior  Neumann  problem.  For  both  the  exterior  Dirichlet  and  Neumann  problems  a 
radiation  boundary  condition  at  00  is  required  to  guarantee  the  uniqueness  of  the  solution. 
This  condition  is  given  by 

lim  r1/2  =  0  .  (4) 

r— »oo  \or  / 

By  developing  local  or  global  boundary  conditions  that  are  applicable  at  a  finite  distance 
R  from  the  scatterer,  it  is  possible  to  replace  (4)  with  such  conditions  and  solve  the 
finite  domain  problem.  Although  the  discrete  problem  has  matrices  which  are  sparse,  the 
condition  number  increases  with  the  inverse  of  the  square  of  the  mesh  size  as  it  approaches 
zero  and  eventually  produces  ill-conditioned  systems  especially  as  k ,  the  wave  number, 
increases.  Instead  we  prefer  an  integral  equation  formulation  which  yields  smaller  matrices 
(unknowns  are  only  on  7  and  not  exterior  to  7),  although  unfortunately  these  matrices  are 
dense. 

We  define  the  field  <f>*0  :  R2  \  {z0}  — »  Cl  of  a  unit  charge  located  at  the  point  ro  €  R2  by 
the  formula 

#w(*)  =  ffo(*ll*-*o||),  (S) 

where  Hq  denotes  the  Hankel  function  of  order  zero.  We  define  the  field  <f>*Q  h  of  a  unity 
dipole  located  at  x0  and  oriented  in  the  direction  h  €  R2  by  the  formula 

#.,»(»)  =  -g»(*ll«  -  «oll)  •  <6) 

For  a  continous  function  a  :  [0,X]  — »  C1,  the  potential  of  a  single  layer  of  density  a  on  a 
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curve  7  is  a  mapping  :  R2  — ►  C1  defined  by  the  formula 


PL(x)  =  f 

Jo 


and  the  potential  of  a  double  layer  of  density  a  on  a  curve  7  is  a  mapping  Pla  :  R2 
defined  by  the  formula 

tL 

ph(x)  =  /  <!>**), 

Jo 


(7) 
C1 

(8) 


Co83  show  that  by  deriving  appropriate  jump  conditions  as  x  — >  7,  the  function  a  satisfies 
the  integral  equation 

2ia(x)  -  Pl,(x)  =  -f(x)  (9) 


Thus,  the  TM  problem  is  solved  by  solving  the  integral  equation  (9),  and  the  resulting 
field  is  generated  by  evaluating  the  integral  (8). 

Similarly,  for  the  TE  problem  (Neumann  case),  using  the  single-layer  potential  (7)  we  must 
solve 

2**(»)  +  ^pk<r(x)  =  9{x)  (10) 

Consider  a  function  VP  satisfying  Helmholtz’s  equation  (1)  find  the  radiation  condition  (4). 
In  polar  coordinates  it  may  be  expanded  as 

»(*)-  f)  0 (11) 

m=- 00 

where  Hm  denotes  the  Hankel  function  of  order  m.  Asymptotic  forms  of  the  radiation  field 
can  be  developed.  First  define 

F*,X0(6)  =  lim  *(<*  +  (12) 

<-►00  y2 
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where  x  =  (cos  9,  sin  0).  Substituting  (11)  into  (12)  and  using  the  asymptotic  decay  of  the 
Hankel  functions,  we  may  write 

00 

JW«)=  £  (13) 

m=— oo 

which  provides  an  explicit  expression  for  FyfXo  via  the  coefficients  and  we  will  refer 
to  Fq/tx0  85  the  asymptotic  representation  of  the  field  $  with  the  origin  at  x0. 

In  order  to  solve  eq.  (9)  or  (10)  numerically  using  Nystrom’s  method,  it  is  necessary 
to  develop  quadrature  formulas  that  can  handle  the  logarithmic  singularity  of  the  Hankel 
function.  Roughly  fourth-order-convergent  quadrature  formulas  are  derived  (Ro85a,  Mu89) 
by  starting  with  the  Euler-Maclaurin  formula  with  the  singular  point  (x  =  0)  removed. 


That  is,  let 


f  f(x)dx  =  h[^2f(xi) 
Jo 


where  h  —  \/n  and  X{  —  i/n.  To  correct  for  the  singular  point,  a  concentration  of  points  of 
the  form  Yfj=i  ^if(Xj)iS  introduced  in  the  first  interval,  where  xj  =jh/6iovj  =  1,2,..., 6. 
The  derivative  term  is  approximated  by  the  one-sided  difference  formula 

/'(*»)  =  -^(3/(x„)  -  4/(x„-!)  +  /(x„-2))  (15) 

Combining  these  terms  yields 

['  f(*)dx  =  h(T  f(Xi)  +  V  Xjfixj)  -  /(x„)/2 

Jo  .=i  i=i 

-  ^(/(x„-2)  -  4/(x„-,)  +  3/(x„)))  (16) 

The  unknown  weights  A j  (j  =  1,2,  ...,6)  are  determined  by  solving  the  6x6  linear  system 
that  results  when  Eq.  (16)  is  assumed  exact  for  the  following  candidate  functions  f(x): 
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1  ,x,x2,logx,xlogx,  and  x2logx ,  using  analytic  integration  rules.  Once  computed,  the 
quadrature  weights  can  be  stored  and  looked  up  numerically  when  needed.  The  proof  of 
approximately  fourth-order  convergence  can  be  found  in  Ro85a. 

Fig.  2  shows  the  approximately  fourth-order  convergence  in  several  electromagnetic  scat¬ 
tering  solutions.  The  relative  errors  are  calculated  by  a  method  which  is  described  in 
Mu89.  The  method  essentially  amounts  to  making  the  code  numerically  verify  Green’s 
second  identity,  and  is  a  very  good  way  to  check  the  accuracy  of  the  computations.  As 
can  be  seen  on  the  figure,  accuracy  at  10  points  per  wavelength  is  at  least  one  percent, 
and  can  be  much  better.  Dramatically  increased  accuracy  can  be  obtained  by  using  more 
points.  This  is  important,  because  integral-equation  methods  commonly  used  today  (most 
of  which  axe  method-of-moments  implementations)  do  not  have  nearly  as  favorable  con¬ 
vergence  properties,  and  some  of  them  even  diverge. 
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G4.  Approaches  to  the  Three-Dimensional  Problem 

The  extension  of  the  FMM  to  three  dimensions  is  non-trivial  because  of  the  difficulty 
of  working  with  spherical  Hankel  functions.  The  essential  idea  is  that  one  replaces  the 
computation  of  the  field  induced  by  a  large  number  of  monopoles  or  dipoles  by  the  com¬ 
putation  of  the  field  induced  by  a  new  computational  element  which  is  a  spherical  Hankel 
function  expansion.  The  cost  of  evaluating  this  expansion  is  much  cheaper  than  the  direct 
evaluation  of  the  field  induced  by  all  the  mono-poles  or  dipoles.  The  cost  of  evaluating  the 
matrix- vector  products  are  reduced  when  the  accuracy  requirements  are  small  because  the 
number  of  terms  in  the  spherical  Hankel  function  expansion  can  be  minimized. 

In  three  dimensions,  the  formulation  in  equations  (1-5)  of  Section  Gl  is  modified  to  use 
the  three  dimensional  free-space  Green’s  function  given  by 

,  exp(tfc|r-r'|) 

G(*|r-r|)  =  -4ir|r_vl 

In  addition,  the  curve  C  in  equation(4)  now  becomes  a  surface  in  three  dimensional  space. 

The  source  distribution  ip( r)  on  C  can  be  approximated  by  linear  functions  over  each 
triangular  element  which  subdivides  the  surface  C,  i.e., 

n 

where  Ln  is  the  hat  function  at  the  point  r„,  the  center  of  the  n-th  triangular  element. 

The  singularities  in  the  integral  in  equation(4)  at  the  point  r  =  r'  may  be  handled  by  the 
quadrature  formulas  in  Mu89. 

Once  equation(4)  is  discretized  using  the  above  expansion  and  quadrature  formulas,  a 
linear  system  results.  The  1-th  component  of  this  system  is  identical  to  the  field  evaluated 
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at  the  point  r i  which  is  induced  by  the  monopoles  or  dipoles  around  the  surface  C.  The 
fast  multipole  method  replaces  this  matrix-vector  operation  by  a  computational  element 
induced  by  these  monopoles  or  dipoles.  Such  a  formula  may  be  derived  by  using  a  Poisson 
formula  for  Helmholtz’s  equation  on  a  sphere  5  given  by 

I(x)  =  J  I(s)H(x,s)ds 


with 

Here  hn^  is  the  first  spherical  Hankel  function,  Pn  is  the  Legendre  polynomial  of  order  n, 
6  is  the  angle  between  Ox  and  Os,  r  =  |x|,  a  =  |s|  and  x  is  a  three  dimensional  point. 
Using  an  appropriate  integration  formula  (e.g.,  Ro85a)  and  truncating  the  above  series  to 
quarantee  a  prescribed  accuracy  requirement  leads  to  the  desired  computational  element 
which  may  be  used  in  the  FMM  to  reduce  the  number  of  calculations. 
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G5  Numerical  treatment  of  corner  singularities. 

The  numerical  methods  that  we  have  employed  to  solve  the  integal  equations  of  scattering 
theory  assume  the  boundaries  of  the  scatterer  are  smooth.  For  trailing  edges  of  airfoils  or 
corner  singularities,  we  must  apply  different  quadrature  formulas.  One  approach  is  given 
below.  Assume  the  parameterization  variable  is  s  with  s  =  0  at  the  comer  and  choose  a 
uniform  mesh,  h.  Then  we  write 


,1  I  N 

/  g(x)  dx  =  V  a'jg(xj)  +  h  V  ff(xj). 
Jo  ;= o 


Suppose  the  corner  singularity  can  be  approximated  by  a  function  like 


g{x)  —  do  +  d\Xa  +  CE2X  + 


where  0  <  cr  <  1  and  a  is  a  function  of  the  comer  angle  or  the  trailing  edge  angle.  The 
unknown  values  aj  in  equation  (1)  can  be  determined  by  assuming  the  formula  is  exact 
for  the  trial  functions  g(x)  =  1,**, *,..., H,  in  addition,  the  function  g(x)  has  a 
logarithmic  singularity  due  to  the  Hankel  function  kernel  in  the  integral  equation,  we  must 
include  in  the  trial  functions  g(x)  those  of  the  form  lnx,xlnx,x2lnx, . . .  as  in  Ro85a, 
Mu89,  and  Mu90.  For  example,  with  a  =  1/2  and  l  =  2,  the  weights  in  equation  (1)  are 
=  .195h,  o.\  =  1.61  h,  and  a'2  =  .695 h.  This  can  be  seen  by  letting  =  hao,  a[  =  hai, 
and  a'2  =  ha 2  +  h/2  and  performing  the  integration  over  the  interval  [0,2hj.  For  g  —  1, 
g  =  x*,  and  g  =  x,  we  obtain 


ha 0  +  hf  x  +  hot 2  =  2 h 


(3a) 


haxh 0  +  ha2(2h)ff  =  (2  h)'+i/(<r  +  1) 


(36) 
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(36) 

Simplifying  equation  (3)  gives 

<*0  +  +  &2  =  2 

(4a) 

+2ffa2  =  2ff+7(a  +  l) 

(46) 

(*1  +  2a2  —  2. 

(4c) 

Solving  equation  (4)  yields 

ao  =  r 

(5a) 

<*i  =  2  —  2r 

(56) 

ot2  —  r 

(5c) 

where 

r  =  [2  —  2<r+1/(<r  +  l)]/(2  —  2ff). 

(5  d) 

Substituting  a  =  1/2  into  equation  (5d)  gives  r  =  0.195,  and  the  result  ioIIows. 
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Fig.  1  Radar  Cross  Section  per  incident  wavelength  (in  dB)  for  first  and  second-order  impedance 
boundary  conditions,  for  a  perfectly  conducting  elliptical  cylinder  (axes  a  —  2,6  =  1; 
kb  =  5)  having  a  thin  coating  of  thickness  6  =  0.005.  The  coating  has  material  parameters 
e  =  n  =  14  + 1.4*.  The  figure  also  shows  a  perfectly  conducting  elliptical  cylinder  without 
a  coating,  for  comparison. 


RELATIVE  ERROR 


NUMBER  OF  POINTS  PER  WAVELENGTH 


Fig.  2  Relative  error  vs.  number  of  points  per  wavelength  on  scatterer  boundary,  for  various 
two-dimensional  conducting  scatterers.  The  reference  line  shows  the  slope  of  fourth-order 
convergence.  The  quadrature  formulas  used  are  approximately  fourth-order  convergent. 
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pole  Method  for  Scattering  from  Electrically  Large  Objects,”  accepted  for  presentation  at 
the  1991  National  Radio  Science  Meeting/  IEEE/AP-S  Symposium,  University  of  Western 
Ontario,  London,  Ontario,  Canada,  June  24-28. 

[2]  Engquist,  B.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “New  Methods  for 
Developing  Higher-Order  Impedance  Boundary  Conditions  on  Curved  Surfaces,”  accepted 
for  presentation  at  the  1991  National  Radio  Science  Meeting/  IEEE/AP-S  Symposium, 
University  of  Western  Ontario,  London,  Ontario,  Canada,  June  24-28. 

[3]  Engheta,  N.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “  What  is  The 
Fast  Multipole  Method  (FMM)  and  How  can  it  Help  You  Solve  Electromagnetic  Scattering 
Problems  More  Effectively?”,  accepted  for  presentation  at  the  Progress  in  Electromagnetics 
Research  Symposium  (PIERS),  Cambridge,  Massachusetts,  July  1-5. 

[4]  Murphy,  W.D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “The  Fast  Multipole  Method 
(FMM)  for  Electromagnetic  Scattering  Problems”,  presented  at  the  URSI/National  Radio 
Science  Meeting,  Jan.  3-5,  Boulder,  Colorado. 

[5]  Murphy,  W.D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “Solving  the  Resonance  Prob¬ 
lem  in  electromagnetic  Scattering  Computation”  ,  presented  at  the  URSI/National  Radio 
Science  Meeting,  Jan.  3-5,  Boulder,  Colorado. 
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[2]  En91a:  Engquist,  B.,  W.  D.  Murphy,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1991),  “Higher 
Order  Impedance  Boundary  Conditions  for  Electromagnetic  Scattering,”  to  be  submitted 
to  Journal  of  Applied  Physics. 

[3]  Mu89:  Murphy,  W.  D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1990),  “Solving  the  Resonance 
Problem  in  Electromagnetic  Scattering  Computation”,  J.  Appl.  Phys.  67  (10),  15  May 
1990,  6061-6065. 

[4]  Mu90:  Murphy,  W.  D.,  V.  Rokhlin,  and  M.  S.  Vassiliou  (1989),  “Numerical  Second- 
Kind- Equation  Solutions  of  Electromagnetic  Scattering  Problems”,  Electronics  Lett.  25, 
643.  (This  is  a  relevant  publication  produced  under  a  related  FY1988  IR&D  project). 
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THE  FAST  MULTIPOLE  METHOD  FOR  SCATTERING 

FROM  ELECTRICALLY  LARGE  OBJECTS 

N.  Engheta*,W.D.  Muiphyt,  V.  Rokhlin*,  and  M.  S. 

Vassiliout 

tDept.  of  Electrical  Engineering,  University  of  Pennsylvania, 
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tRockwell  International  Science  Center,  1049  Camino  Dos 

Rios,  Thousand  Oaks,  CA  91360 

*Dept.  of  Computer  Science,  Yale  University,  New  Haven  CT 

06520 

The  Fast  Multipole  Method  (FMM)  was  developed  by  Rokhlin  (7. 
Comp.  Phys.  60,  187-207,  1985;  Yale  Univ.  Research  Repor' 
YALEU/DCS/RR-440, 1985)  to  solve  acoustic  scattering  problems  very 
efficiently.  We  have  modified  and  adapted  it  to  electromagnetic  scattering 
problems  in  two  dimensions. 

Consider  a  2-D  closed  conducting  scatterer  with  n  nodes  on  its 
boundary.  Divide  the  boundary  into  p  segments,  each  containing  nip 
nodes.  Instead  of  calculating  n2  interactions  among  n  current  sources  on  the 
boundary,  consider  each  segment  to  be  a  cluster  of  nip  sources.  For 
segments  that  arc  close  together,  the  exact  interactions  must  be  calculated. 
For  segments  sufficiently  far  apart,  however,  we  may  do  the  following:  we 
combine  the  sources  in  each  segment,  approximating  their  radiation  fields 
by  the  first  N  multipoles.  We  describe  each  segment  via  an  equivalent 
source  located  at  the  segment's  center.  We  can  calculate  the  contribution  of 
each  such  equivalent  source  to  the  field  at  the  center  of  any  sufficiently 
separated  receiver  segment,  and  in  that  receiver  segment  we  can  use  a 
Taylor  expansion  to  obtain  the  field  at  all  the  individual  nip  nodes.  The 
radiation  field  at  any  particular  node  on  the  scanerer  boundary  is  the  sum  of 
contributions  of  N  multipoles  of  each  of  the  far-away  segments,  and  ‘he 
direct  contribution  of  very  close  segments. 

The  method  reduces  the  operation  count  for  solving  the 
Magnetic-Field  Integral  Equation  (MFIE)  from  0(n3)  for  Gaussian 
elimination  to  0(n4/3)  per  conjugate-gradient  iteration.  It  has  proven  useful 
in  calculating  the  scattering  from  electrically  large  objects  difficult  to 
compute  by  many  other  methods. 


'This  work  was  supported  by  Air  Force  Office  of  Scientific  Research 
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NEW  METHODS  FOR  DEVELOPING  HIGHER  ORDER 
IMPEDANCE  BOUNDARY  CONDITIONS  ON  CURVED 
SURFACES1 

B.  Engquist^W.D.  Murphyt,  V.  Rokhlin*,  and  M.  S. 

Vassiliout 

tDept  of  Mathematics,  University  of  California  at  Los  Angeles, 
Los  Angeles,  CA  90024 

tRockwell  International  Science  Center,  1049  Camino  Dos 
Rios,  Thousand  Oaks,  CA  91360 

♦Dept,  of  Computer  Science,  Yale  University,  New  Haven  CT 
06520 

We  have  developed  techniques  for  computing  electromagnetic 
scattering  from  closed  2-D  conductors  coated  with  multiple  thin  layers  of 
possibly  lossy  dielectric  and/or  magnetic  material.  Using  Fourier  integral 
techniques,  we  derive  higher-order  impedance  boundary  conditions  of 

0(82)  and  0(84)  in  the  thickness  parameter  8.  We  develop  second-kind, 
weakly  singular  integral  equations.  These  integral  equations  are  discretized 
and  solved  using  Nystrttm’s  method  and  approximately 
fourth-order-convergent  quadrature  formulas.  Solutions  compare 
favor'ably  with  analytical  results.  We  have  used  our  techniques  to  study  the 
effects  of  layer  thickness,  body  geometry,  and  incidence  angle  on  the 
scattered  fields. 

mis  work  was  partially  supported  by  Air  Force  Office  of  Scientific 
Research  Contract  Number  F49620-89-C-0048 
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We  describe  the  Fast  Multipole  Method  (FMM)  developed  by  Rokhlin 
(J.  Comp.  Phys.  60,  187-207,  1985;  Yale  Univ.  Research  Report 
YALEU/DCS/RR-440, 1985)  to  reduce  the  computation  time  required  for 
solving  scattering  problems  modeled  in  the  form  of  integral  equations.  The 
method  was  developed  first  for  acoustics.  We  have  modified  and  adapted  it 
to  electromagnetic  scattering  problems  in  two  dimensions.  Briefly,  the 
method  relies  on  dividing  the  boundary  of  the  scatterer  into  segments,  each 
containing  several  nodes,  and  performing  exact  computations  only  for  the 
interactions  of  neighboring  segments.  When  segments  are  sufficiently  far 
from  each  other,  approximations  are  used,  expanding  “source  segments” 
into  multipoles. 

Our  present  computer  code  reduces  the  operation  count  for  solving  the 
Magnetic-Field  Integral  Equation  (MFIE)  from  0(n3)  to  0(n4/3)  per 
conjugate-gradient  iteration.  It  has  achieved  dramatic  results  in  solving 
scattering  from  electrically  large  objects.  Hie  implementation  is  presently 
limited  to  computing  scattering  from  closed  conducting  objects  in  two 
dimensions.  Impedance  boundary  conditions  which  have  been  developed 
and  tested  in  a  simpler  computer  code  will  soon  be  merged  into  the  FMM 
implementation  to  allow  computation  of  conductors  coated  with  dielectrics. 
In  addition,  work  is  proceeding  on  a  three-dimensional  implementation. 


iThis  work  was  supported  by  Air  Force  Office  of  Scientific  Research 
Contract  Number  F49620-89-C-0048 
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The  Fast  Multipole  Method  (FMM)  was  developed  by  Rokhlin  (/. 
Comp.  Phys.  60,  187-207,  1985;  Yale  Univ.  Research  Report 
YALEU/DCS/RR-440, 1985)  to  reduce  the  computation  time  required  for 
solving  scattering  problems  modeled  in  the  form  of  integral  equations. 
When  such  integral  equations  are  discretized,  a  full  dense  linear  system  of 
size  nxn  results.  If  direct  methods  such  as  Gaussian  Elimination  are  used 
to  solve  the  system,  0(n3)  arithmetic  operations  are  needed.  If  indirect 
methods  such  as  Generalized  Conjugate  Residual  (GCR)  are  employed, 
0(n2)  arithmetic  operations  are  needed  for  each  iteration.  Rokhlin  has 
shown  that  the  operation  count  may  be  reduced  to  0(n  logn  )per  iteration, 
by  using  the  structure  of  the  Green's  function  and  the  fact  that  only  an 
approximation  to  the  matrix  product  is  needed.  The  algorithm  is  analogous 
to  the  evaluation  of  the  field  created  by  a  charge  and  dipole  distribution 
(hence  the  name  "Fast  Multipole  Method").  The  actual  implementation  thus 
far  uses  a  simpler  algorithm  which  achieves  0(n4/3),  although  0(n  log n )  is 
achievable  theoretically.  For  moderately  large  problems  (n  *  20000),  the 
two  algorithms  yield  similar  CPU  times,  although  for  very  large  problems, 
implementing  the  full  O (n  logn )  algorithm  is  desirable. 

Originally,  Rokhlin’s  work  was  for  acoustic  scattering  problems 
where  a  fluid  scatterer  is  embedded  in  a  two-dimensional  fluid  space.  This 
mathematical  problem  is  formulated  as  a  coupled  system  of  integral 
equations  between  the  interior  and  exterior  problems.  Rokhlin  extended 
these  ideas  also  to  two  and  three-dimensional  potential  theory  (Laplace's 
equation)  where  even  faster  computations  are  possible  (0(n)  operations). 

We  review  some  of  these  results,  and  extend  their  application  to  the 
area  of  electromagnetic  scattering.  Initial  results  with  the  0(n4/3)  algorithm 
have  been  obtained  for  TM  scattering  from  arbitrarily  shaped 
two-dimensional  closed  conductors.  The  FMM,  combined  with 
fourth-order  convergent  quadrature  formulas  in  Nystrom's  method 
(Murphy,  Rokhlin  and  Vassiliou,  Electronics  Letters  25,  543-644,  1989) 
represents  a  very  fast  and  accurate  method  to  solve  electromagnetic 
scattering  problems  from  electrically  large  objects. 
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The  "resonance  problem"  considered  here  is  that  at  certain  values  of 
the  wavenumber  k ,  the  second-kind  integral  equation  for  solving  scattering 
problems  can  become  extremely  ill-conditioned.  This  affects  both  the 
accuracy  and  speed  of  numerical  solutions.  The  convergence  of  iterative 
methods  in  particular  is  adversely  affected  by  large  condition  numbers.  We 
consider  the  second-kind  integral  equation,  derived  using  double-layer 
potentials,  for  TM  scattering  from  a  conductor.  We  solve  an  exterior 
Dirichlet  problem.  It  has  nonunique  solutions  for  values  of  k  at  which  the 
corresponding  homogeneous  interior  Neumann  problem  has  a  nontrivial 
solution.  These  values  of  k  are  the  resonant  wavenumbers. 

Numerically,  we  discretize  the  integral  directly  using  quadrature 
formulas  (Nystrom's  method).  We  employ  fourth-order  convergent 
quadrature  formulas  which  handle  the  logarithmic  singularities  in  the 
problem  (Murphy,  Rokhlin,  and  Vassiliou,  Electron  Lett.  25, 643,  1989). 
We  tested  the  solution  procedure  at  resonant  k' s  for  circular  and  elliptical 
scatterers  (roots  of  derivatives,  respectively,  of  Bessel  functions  and 
modified  Mathieu  functions).  We  found  very  large  condition  numbers  for 
the  discrete  matrices  (up  to  O(107)),  generally  leading  to  poor  solutions. 

We  apply  two  approaches  to  alleviate  the  resonance  problem.  The 
first  is  to  use  a  different  integral  equation,  based  on  both  single  and 
double-layer  potentials,  analogous  to  the  combined-field  equation  (CFEE). 
This  leads  to  low  condition  numbers,  and  good  solutions,  at  resonant  k.. 
The  second  method  is  to  use  the  original  second-kind  integral  equation, 
introduce  a  small  imaginary  part  in  k,  and  extrapolate  back  to  the  real  axis. 
Solutions  obtained  by  the  two  methods  are  in  excellent  agreement. 

By  solving  the  resonance  problem,  we  ensure  that  fast  and  accurate 
solutions  are  obtainable  at  any  arbitrary  wavenumber. 

1 77215  work  was  partially  supported  by  Air  Force  Office  of  Scientific 
Research  Contract  Number  F49620-89-C-0048 
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ABSTRACT 


The  Fast  Multipole  Method  (FMM)  was  developed  by  Rokhlin  to  solve  acoustic  scat¬ 
tering  problems  very  efficiently.  We  have  modified  and  adapted  it  to  the  second-kind- 
integral-equation  formulation  of  electromagnetic  scattering  problems  in  two  dimensions. 
The  present  implementation  treats  the  Dirichlet  (TM)  problem  for  two-dimensional  closed 
conducting  objects  of  arbitrary  geometry.  Consider  a  scatterer  with  n  nodes  on  its  bound¬ 
ary.  Divide  the  boundary  into  p  segments,  each  containing  n/p  nodes.  Instead  of  calcu¬ 
lating  n2  interactions  among  n  current  sources  on  the  boundary,  consider  each  segment  to 
be  a  cluster  of  n/p  sources.  For  segments  that  are  close  together,  the  exact  interactions 
must  be  calculated.  For  segments  sufficiently  far  apart,  however,  we  may  combine  the 
sources  in  each  segment,  approximating  their  radiation  fields  by  the  first  N  multipoles. 
We  describe  each  segment  via  an  equivalent  multipole  located  at  the  segment’s  center.  We 
can  calculate  the  contribution  of  each  such  equivalent  expansion  to  the  field  at  the  center 
of  any  sufficiently  separated  receiver  segment,  and  in  that  receiver  segment  we  can  use  a 
partial  wave  expansion  to  obtain  the  field  at  all  the  individual  n/p  nodes.  The  radiation 
field  at  any  particular  node  on  the  scatterer  boundary  is  the  sum  of  contributions  of  N 
multipole  expansions  of  each  of  the  far-away  segments,  and  the  direct  contribution  of  very 
close  segments.  The  method  reduces  the  operation  count  for  solving  the  Magnetic-Field 
Integral  Equation  (MFIE)  from  0(n3)  for  Gaussian  elimination  to  0(nA/:i)  per  conjugate- 
gradient  iteration.  We  also  present  a  simple  technique  for  accelerating  convergence  of  the 
iterative  method:  “complexifying”  k,  the  wave  number.  This  has  the  effect  of  bounding  the 
condition  number  of  the  discrete  system;  consequently,  the  operation  count  of  the  entire 
FMM  (all  iterations)  becomes  0(n4/3).  We  present  computational  results  for  moderate 
values  of  ka ,  where  a  is  the  characteristic  size  of  the  scatterer. 
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I.  INTRODUCTION 


Electromagnetic  scattering  comprises  an  important  class  of  problems  in  physics  and 
engineering.  It  is  desirable  to  have  efficient  techniques  to  compute  scattering  accurately.  In 
this  paper  we  present  the  first  application  of  Rokhlin’s  Fast  Multipole  Method  (FMM)  to 
the  problem  of  electromagnetic  scattering  from  two-dimensional  closed  conducting  objects 
of  arbitrary  geometry. 

In  many  two-dimensional  scattering  problems,  it  is  customary  to  reduce  the  scalar 
Helmholtz  equation  to  a  second-kind  integral  equation.  The  resulting  integral  equation 
can  generally  be  treated  using  various  numerical  techniques1,2.  One  of  the  standard  meth¬ 
ods  for  this  numerical  treatment  of  scattering  problems  is  to  discretize  the  second  kind 
integral  equation  using  an  appropriate  quadrature  formula  (Nystrom’s  method)1,3.  Such 
discretization  leads  to  systems  of  linear  algebraic  equations  which  may  be  solved  by  Gaus¬ 
sian  elimination  or  iterative  methods  such  as  conjugate  gradient  or  generalized  conjugate 
residual  (GCR)4.  These  iterative  methods  require  the  full  dense  matrix  to  operate  on 
a  sequence  of  recursively  generated  vectors.  Consequently,  the  operation  count  is  0{n2) 
where  n  is  the  dimension  of  the  matrix.  There  have  been  many  successful  efforts  to  reduce 
the  operation  count  and  storage  requirements  and  to  introduce  “fast’’  algorithms5, 6,1,8 . 
The  FMM  is  a  particularly  promising  method  among  these.  In  the  FMM,  the  operation 
count  for  each  iteration  is  reduced  to  0(?i4/,;5),  which  is  significantly  smaller  than  0(n2), 
especially  for  large  ?z(>  10,000).  This  algorithm  can  be  further  improved  to  one  that  has 
an  operation  count  of  0{n  log  n)  per  iteration.  However,  we  have  not  yet  implemented  the 
fastest  algorithm.  When  these  algorithms  are  combined  with  a  GCR  or  conjugate  gradient 
algorithm,  the  resulting  procedure  only  requires  a  small  number  of  iterations  to  converge 
to  a  solution  to  the  scattering  problem.  This  is  the  case  even  at  resonance  frequencies,  if 
the  method  of  “complexification”  is  applied  (see  discussion  below). 

The  purpose  of  the  present  paper  is  to  explain  this  algorithm  in  simple  terms  and 
to  explore  its  application  to  electromagnetic  scattering.  Rokhlin’s  FMM,  which  was  first 
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employed  for  Poisson’s  equation5  and  acoustic  scattering6  in  two  dimensions  has  been 
extended  to  three  dimension  for  Poisson’s  equation9  and  is  currently  being  extended  for 
Helmholtz’s  equation. 
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II.  PROBLEM  STATEMENT 


Consider  a  two-dimensional  conducting  body  whose  axis  is  aligned  with  the  z  co¬ 
ordinate  axis.  A  monochromatic  electromagnetic  wave  incident  on  this  structure  with  an 
electric  field  vector  parallel  to  the  axis  of  the  body  is  referred  to  as  the  transverse-magnetic 
(TM)  case.  The  incident  and  scattered  fields  both  satisfy  the  following  Helmholtz  equation 

A  Ez  +  k2Et  =  0.  (1) 

The  boundary  condition  for  equation  (1)  is  that  the  total  E  field  vanish  on  the  surface  of 
this  conductor,  i.e., 

E\ot  -  0  on  C  (2) 

or  more  explicitly 

E\nc  +  E?cat  =  0  on  C.  (3) 

The  above  Helmholtz  equation  can  be  reduced  to  the  following  second-kind  integral  equa¬ 
tion 

1'(r)  +  2^9?.W^;)r'l)i’(r')d1'  =  -2E‘"c(r)  (4a) 

£r'(r>= (4,,) 

where  r  and  r'  are  both  on  the  boundary  C,  and  G  is  the  free-space  Green’s  function  in 
two  dimensions,  i.e., 

G(fc|r  -  r'|)  =  iH<‘>(k|r  -  r')/4  (5) 

with  -  r'|)  being  the  Hankel  faction  of  the  first  kind  of  order  zero.  See  Hef. [10} 

for  a  derivation  of  equation  (4).  Equation  (4),  the  TM  case,  is  often  referred  to  as  the 
Dirichlet  problem  in  the  mathematical  literature1,10.  If  we  discretize  the  boundary  into  n 
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points,  then  the  above  integral  equation  (4a)  is  converted  to  the  following  linear  system 
(via  Nystrom’s  method1 ): 


*M+2£au¥(.5)  =  — 2ET(r,)  (6) 

j  =  l 

where  the  matrix  A  =  (AtJ)  is  n  x  n  and  the  vectors  (^(rj))  and  (E^"c(rj))  are  col¬ 
umn  vectors  having  n  rows.  Applying  normal  matrix  multiplication,  A$  requires  0(n2) 
operations.  The  FMM  algorithm  reduces  this  to  0(n4/3)  or  ultimately  to  O(nlogn). 
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III.  RAPID  SOLUTION  OF  INTEGRAL  EQUATIONS 


The  detailed  mathematics  behind  the  FMM  is  presented  in  [5]  and  [6].  The  develop¬ 
ment  is  quite  complex.  Below,  we  offer  a  simplified  version,  with  more  physical  intuition 
relevant  to  electromagnetic  scattering. 

Consider  n  nodes  on  the  boundary  of  the  scatterer.  Divide  the  boundary  into  p  equal 
segments,  where  2  <  p  <  n.  In  each  segment,  there  are  njp  nodes.  If  the  length  of  the 
boundary  is  L,  each  segment  has  length  L/p.  The  center  of  each  segment  is  located  at 
Z{(i  =  1,2,..., p).  In  scattering  problems,  each  node  can  be  treated  conceptually  as  if  it. 
were  a  source  of  radiation. 

If  we  have  sources  within  a  finite  region  of  space,  the  radiation  emitted  from  these 
sources  in  the  far  zone  can  be  approximated  using  a  collection  of  multipoles  located  at  the 
center  of  the  region5,6,11 .  The  multipole  approximation  converges  rapidly  outside  any  circle 
D  containing  all  sources  and  separated  from  D  by  at  least  one  wavelength.  In  fact,  once 
a  sufficiently  large  number  of  multipoles  is  included,  the  accuracy  of  the  approximation 
increases  superalgebraically  (faster  than  any  negative  power  of  N)6. 

Consider  each  segment  on  the  boundary  as  a  cluster  of  n/p  sources.  The  sources  in  each 
segment  are  treated  as  a  single  aggregate  source,  and  the  radiation  field  of  that  equivalent 
source  is  approximated  using  the  first  N  multipoles  located  at  the  center  of  the  segment. 
For  each  pair  of  sufficiently-separated  segments,  the  radiation  of  the  N  multipoles  of  one 
segment  can  be  represented  as  an  analytical  partial  field  expansion  around  the  center  of  the 
other  segment.  Then  from  this  information,  the  field  at  the  other  nodes  of  that  segment 
can  be  evaluated  using  equation  (16).  For  nearby  segments,  the  direct  contribution  must 
be  calculated  to  evaluate  the  radiation  field.  The  radiation  field  at  any  particular  node 
on  the  boundary  is  the  sum  of  the  contributions  of  N  multipoles  of  each  of  the  far-away 
segments  and  the  direct  contribution  of  the  nearby  segments.  Reference  [C]  considers  the 
precise  mathematical  description  of  the  process. 
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To  illustrate  the  above  verbal  description  mathematically,  let  us  consider  the  field  of  a 
two-dimensional  magnetic  surface  current  distributed  over  a  two-dimensional  body.  Since 
we  have  an  exterior  Dirichlet  problem,  a  double-layer  potential  is  needed.  For  the  two- 
dimensional  TM  problem,  a  magnetic  current  density  at  any  point  of  the  boundary  is  a 
vector  in  the  x-y  plane  tangent  to  the  body.  The  electric  field  of  such  a  current  distribution 
is  given  by 

B ■“*(«  =  V  x  G(k|/T—  ^l)K(p') dl*  (7) 

where  K(p')  is  the  magnetic  surface  current  density,  p  =  (p,  #),  and  C  is  the  contour  of 
integration.  In  two  dimensions  equation(7)  can  be  written  as 


E 


scat 

z 


dG(k\p-p'\) 

dn(p') 


|K(p')|dl' 


(8) 


which  is  identical  with  equation(4b)  if  |K(p')|  is  taken  to  be  'i’(p').  Substituting  equa¬ 
tion^)  into  (8),  we  get 

Ef-M  =  (i/4)  Jc  (9) 

The  Hankel  fuction  H^^p  —  p'\)  can  be  expanded  in  terms  of  higher  order  Hankel 
functions: 

OO 

^,)W-/r'l)=  E  B^(b)Jm(b')ext (im(e-0'»  (10) 

m=  — oo 

Substituting  equation  (10)  into  equation  (9)  yields 


o° 

Elcat(p,0)=  H^(kp)exp(imd)  /  (i/4) 

m  ——oo  JC 


dJm(kp')exp(—im6') 

dn(p',6') 


\K(p',0')\d\'  (11) 


This  can  be  regarded  as  the  multipole  expansion  of  the  source  K (p',61).  For  a  discretized 
source  at  n  points  located  at  x'}  =  ( p) ,  8'} )  (j  =  1, 2, . . .  ,n)  over  a  segment  of  the  boundary, 
equation(ll)  reduces  to 


E 


scat 


.  A  v~v  ,  dJm(kp'i)exp(—im6'i)  ... 

(/>.*)=  E  D1/4) -  "dn(x/) - (12) 


m  =  —  oo  j—1 
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where  A/'-  is  the  discretized  element  of  arc  length  containing  the  source  |K(xj)|.  For  a 
given  accuracy,  we  can  truncate  the  infinite  sum  in  equation(12)  at  Ar,  and  thus  calculate 
the  first  N  multipoles  of  the  source. 

Historically,  the  FMM  algorithm  was  first  applied  to  Poisson’s  equation  for  n  point 
charge  lines  at  locations  x,(i  =  1,2, ...,n)  with  strengths  k;.  This  is  mathematically 
equivalent  to  solving  the  equation 

n 

A<f>  =  ^  6(x  —  Xi)ni  (13) 

i=i 

where  6(x)  is  the  Dirac  delta  function  and  x  and  x,  are  two-dimensional  points.  The 
solution  to  equation  (13)  is 

n 

^(x) =  53 Ki  los^x  ~  x»i)/(27r)  (14) 

t=i 

If  we  evaluate  equation  (14)  at  each  point  x,(?  =  1,2,  ..,.,n),  then  this  computation  re¬ 
quires  0(n2)  operations.  However,  if  large  numbers  of  particles  are  combined  into  single 
computational  elements,  then  this  operation  »  xint  can  be  reduced  if  an  approximate  an¬ 
swer  (to  a  specified  accuracy)  is  desired.  When  a  cluster  of  particles  is  “far  away”  from  a 
particular  point,  then  the  potential  of  the  cluster  is  approximated  by  the  potential  induced 
by  a  single  computational  element  located  inside  the  cluster12.  In  the  FMM  algorithm  the 
computational  element  is  a  Laurent  expansion  centered  at  a  circle  containing  the  cluster 
of  particles.  Given  a  cluster  of  charges  located  at  points  z,(i  =  1,2, . . .  ,«c),  the  expansion 
is  given  by 

nc  p 

<j>(z)  =  Re(^2 l°g(*  -  zi))  ~  i?e(a0log(2  -  z0)  +  53°*/^  “  5o)*)  (15) 

j=l  k—l 

Here  p  is  the  order  of  the  multipole  and  the  ajts  are  coefficients  chosen  so  that  the  truncated 
series  is  an  accurate  approximation  of  the  potential.  For  equation  (15),  the  computational 
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effort  is  only  0(p)  operations  and  is  much  less  than  0(nc)  for  the  direct  approach.  The 
region  must  now  be  organized  into  well  separated  points  and  very-near  points.  For  near 
points,  the  direct  evaluation  of  (14)  is  used.  See  Refs. [5,6,12]  for  details  of  decomposing 
the  regions  into  boxes  of  different  sizes.  When  applying  the  FMM  to  Helmholtz’s  equa¬ 
tion  instead  of  Poisson’s  equation  the  expansion  (15)  is  replaced  by  the  standard  Hankel 
function  expansion,  which  is  then  truncated  to  obtain  a  given  accuracy  requirement.  That 
is, 

Eacat(o  e\  =  {  Sm=-ooam^m)(^p)exp(trn0),  if  P  >  a  (16a) 

'  1  P™Jm(kp)  exp(im8),  if  p<a  (16b) 

depending  on  whether  the  calculations  arc  to  be  done  outside  or  inside  the  circle  of  adius 
a6.  Obviously,  a iij  prescribed  accuracy  in  the  series  can  be  guaranteed  by  taking  more 
terms  in  the  expansion  at  the  expense  of  more  CPU  time. 
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IV.  OPERATION  COUNT 


We  illustrate  the  computational  work  required  for  the  FMM  algorithm  using  a  simple 
example.  Consider  p  segments  around  the  boundary  of  the  scatterer.  Assume  that  for 
each  segment,  the  two  adjacent  segments  are  “nearby,”  and  thus  require  direct  calculation 
of  the  radiation  field.  All  other  segments  are  considered  “far  away”  for  this  example  and 
the  multipole  expansion  will  be  used.  The  following  steps  are  taken: 

Step  1:  Find  the  first  N  multipoles  of  sources  in  each  segment.  Since  to  evaluate  each 
multipole,  all  the  sources  are  involved,  the  first  N  multipoles  of  sources  in  each  segment 
require  Nri/p  operations.  For  p  segments,  we  have 

N(n/p)p  =  Nn  (17) 

As  we  mentioned  earlier,  the  number  of  multipoles  used  is  a  function  of  the  accuracy 
needed  in  the  calculation.  This  number  is  typically  proportional  to  kd  where  k  =  2n/X 
and  d  is  the  length  of  the  source  region.  Thus,  N  «  kL/p.  Substituting  this  value  of  N 
in  o  equation  (17)  gives 

Nn  «  k(L/p)n  (18) 


as  the  operation  count  for  this  step. 

Step  2:  For  each  pair  of  “far  away”  segments,  evaluate  the  radiation  fields  of  N  multipoles 
of  one  segment  at  the  center  of  t lie  other.  This  is  an  0(N)  operation  for  each  pair.  The 
number  of  pairs  is  almost  p2.  In  actuality,  for  each  segment,  the  number  of  far  away 
segments  is  p  —  3.  Therefore,  the  number  of  far  away  pairs  is  p(p  —  3),  which  for  large  p  is 
almost  p2.  Therefore,  the  operation  count  for  this  step  is 

Np 2  «  k(L/p)p2  ss  kLp  (19) 
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Step  3:  In  this  step,  for  each  segment,  add  the  contribution  of  N  multipoles  of  any  one 
of  the  far  away  segments  evaluated  at  the  center  of  the  chosen  segment.  For  any  chosen 
segment,  the  number  of  far  away  segments  is  p  —  3,  or  approximately  p,  for  large  p.  Thus, 
this  step  requ:res 

Np  «  k(L/p)p  ss  kL  (20) 


operations. 

Step  4:  Here,  the  radiation  field  is  known  at  the  center  of  each  segment.  The  field  at  the 
other  nodes  in  the  segment  can  be  evaluated  using  a  partial  field  expansion.  For  each 
neighboring  node,  this  is  an  0(N )  operation.  Thus,  for  n/p  —  1  nodes  in  each  segment, 
the  number  of  operations  is 

N(n/p  -  1)  «  Nn/p  (21) 


For  p  —  3  segments,  we  have 


(Nn/p)(p  -  3)  S3  Nn  R3  k(L/p)ri  (22) 

Step  5:  Finally  for  the  nearby  (neighboring)  segments,  the  direct  contributions  must  be 
evaluated.  For  q  sources,  the  number  of  operations  is  q2.  Here,  in  each  segment  there  are 
n/p  sources.  For  a  particular  segment  in  question  and  its  two  near-neighbors,  the  work  is 

(3 n/p)2  ?3  n2/p2  (23) 


For  p  segments,  the  count  is 

(n2/p2)p~n2/p  (24) 

Adding  the  above  five  steps  and  optimizing  with  respect  to  p  the  resulting  expression 
for  total  operation  count  (as  in  Ref.[6]),  the  optimal  count  T  is 

T  =  0(n3/2)  (25) 
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The  operation  count  can  be  further  reduced  by  applying  the  above  procedure  recursively, 
with  each  of  the  nearby  segments  subdivided  with  appropriately  chosen  p'.  The  new 
estimate  so  obtained  shows  that  the  FMM  algorithm  is  0(n4/3).  Our  current  FMM  code 
implements  the  0(n4//3)  algorithm.  In  Ref. [5]  the  above  subdivision  is  used  recursively  to 
obtain  an  0(n)  algorithm  when  applied  to  Poisson’s  equation  (13).  By  reproducing  the 
construction  of  Section  VII  in  Ref.[5]  for  Helmholtz’s  equation  an  0(n  log(n))  algorithm  is 
theoretically  achievable. 

The  essential  feature  of  the  FMM  algorithm  is  that  it  performs  the  matrix-vector 
operation 

A$(m)  (m  =  0,1,2,...)  (26) 

in  0(n4/3)  operations.  Here  the  superscript  m  is  the  iteration  counter.  Note  that  A  is  never 
computed,  so  that  the  algorithm  only  requires  vector  storage  !  The  storage  requirement  is 
0(n4/3),  like  the  operation  count6.  The  actual  solution  to  the  system  (6)  is  obtained  using 
a  conjugate  gradient  or  GCR4  algorithm  in  which  the  most  computationally  extensive  step 
is  that  of  forming  the  vector  expression  (26)  by  the  FMM  technique.  GCR  must  be  used 
on  the  matrix  I  +  2A  or  conjugate  gradient  on  the  normal  equation  because  I  +  2.4  is 
nonsymmetric. 

The  FMM  algorithm  is  in  no  way  restricted  to  the  TM  case  (4).  In  fact,  an  almost 
identical  algorithm  can  be  applied  to  the  TE  case  (Neumann  problem  in  the  mathematical 
literature1,10)  or  to  the  combined  field  integral  equation  (CFIE),  which  does  not  have  a 
resonance  problem  (see,  e.g.,  [3]).  Rather  than  solve  the  CFIE,  we  use  a  “lazy  man’s” 
technique  that  employs  the  method  of  “complexification”  on  equation  (4)  directly.  This 
“trick”,  the  reader  will  soon  see,  works  remarkably  well  and  allows  us  to  avoid  dealing  with 
the  more  complicated  CFIE.  Although  “complexification”  has  been  used  in  electromagnet¬ 
ics  to  satisfy  the  radiation  condition,  here  it  is  used  to  reduce  the  condition  number  of  the 
matrix  A  and  thereby  reduce  the  number  of  iterations  for  convergence  of  the  conjugate 
gradient  algorithm. 
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V.  COMPUTATIONAL  RESULTS 


In  Fig.  la,  we  show  the  comparison  of  the  FMM  algorithm  with  the  analytic  solution 
for  the  scattering  of  a  plane  wave  incident  on  a  circular  cylinder  for  the  case  ka  =  80. 
We  present  our  computational  results  in  a  form  often  used  in  electrical  engineering:  that 
is,  as  plots  of  differential  scattering  cross  section  or  Radar  Cross  Section  (RCS)  versus 
observation  angle,  for  a  given  angle  of  incidence.  The  RCS  is  related  to  the  magnitude  of 
the  electric  field  in  the  far  zone.  In  two  dimensions,  “RCS”  is  something  of  a  misnomer:  the 
more  proper  term  is  “scattering  width,”  but  the  label  RCS  is  commonly  applied  anyway. 
The  two-dimensional  definition  is13 


RCS  =  a  —  2ir  lim  r 

r—+  oo 


(27) 


where  Eszcat  is  the  scattered  field  and  E‘nc  is  the  incident  field.  The  quantity  we  plot  is 
the  ratio  of  a  to  the  wavelength  of  the  incident  wave,  expressed  in  decibels  (dB). 

Fig.  la  shows  the  agreement  of  the  RCS  between  the  two  solutions  is  exact  for  the 
observation  angle  <j>  between  0°  and  40°,  where  most  of  the  rapid  changes  occur.  In  Fig. 
lb,  we  plot  the  entire  FMM  solution  from  0 0  <  <f>  <  180°,  which  also  agrees  exactly  with 
the  series  solution,  but  we  leave  out  the  latter  to  avoid  too  much  congestion  in  the  figure. 
In  this  example,  we  have  used  10  points  per  wavelength  on  the  scattcrer  boundary,  so 
n  =  800.  In  Table  1,  we  list  some  results  from  the  FMM  code  for  various  values  of  k  and  n 
for  scattering  from  two-dimensional  conducting  circular  and  elliptical  cylinders  (a=semi- 
major  axis,  b=semi-minor  axis).  The  expansions  (16)  were  truncated  to  give  an  accuracy 
of  10~4  and  the  convergence  tolerance  for  the  conjugate  gradient  algorithm  was  set  to 
10-3.  In  the  table,  Ntter  denotes  the  actual  number  of  iterations  in  the  conjugate  gradient 
algorithm  to  achieve  a  relative  error  tolerance  of  10-3,  ERR  is  the  actual  relative  error 
tolerance  between  the  last  two  iterations,  and  CPU  is  the  average  CPU  time  in  seconds  on 
a  VAX  6410  computer  for  one  iteration  of  the  conjugate  gradient  FMM  algorithm  using 
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double  precision  arithmetic.  In  most  cases,  for  “complexified”  k-values  an  eror  tolerance 
of  10“4  instead  of  10-3  would  only  add  one  or  two  more  iterations  to  Ar,(er.  However,  for 
“non-complexified”  k-values,  an  error  tolerance  of  10“4  would  require  many  more  iterations 
because  of  the  large  condition  number  involved  in  those  cases.  Since  the  matrix  A  is  never 
explicitly  calculated  in  the  FMM  algorithm,  the  total  CPU  time  for  one  incident  angle  is 
roughly  equal  to  the  product  of  ITER  and  CPU.  To  calculate  the  number  of  points  N\  per 
incident  wavelength  in  Table  1,  we  have  used  an  approximate  formula  for  the  perimeter  of 
an  ellipse,  obtaining 


Note  that  for  some  of  the  cases  in  Table  1,  we  have  used  values  of  k,  the  wave  num¬ 
ber,  that  have  imaginary  parts.  This  was  done  to  test  the  method  of  “complexification”. 
Equation  (4)  has  resonances10  only  at  real  distinct  values  of  k  and  as  k  gets  larger,  root 
clusters  become  more  dense;  at  these  resonance  frequencies,  the  condition  number  of  the 
matrix  A  becomes  large  and  iterative  methods  require  many  more  iterations  to  converge 
to  a  given  tolerance.  We  have  shown  in  Ref.[3]  that  by  moving  k  slightly  into  the  com¬ 
plex  plane(  “complexifying”  k)  the  condition  number  can  be  reduced  by  four  or  five  orders 
of  magnitude  and  consequently,  the  convergence  rate  of  most  iterative  methods  can  be 
greatly  improved.  This  is  clearly  demonstrated  in  Table  1.  Furthermore,  by  using  two 
complex  values  of  k  (say  k  =  100  +  .li  and  k  =  100  +  .2?)  we  have  shown  in  Ref. [3]  that 
the  extrapolated  values  of  the  RCS  to  the  real  axis  are  accurate  to  a  value  in  dB  roughly 
equal  to  the  amount  of  movement  into  the  complex  plane  (in  this  case  .ldB).  Of  course, 
parabolic  extrapolation  would  be  even  more  accurate,  but  would  require  computations  for 
three  complex  values  of  k.  However,  this  is  not  a  severe  problem  when  using  iterative 
methods  because  the  iteration  scheme  for  the  second  value  of  k  can  be  started  with  the 
final  solution  from  the  first  value  of  k  and  so  on.  Thus  two  solutions  for  two  close  values 
of  k  require  only  about  1.1  times  as  much  CPU  time  as  does  one  case.  Complexification  is 
hence  a  viable  and  realtively  cheap  method  for  accelerating  convergence,  in  such  problems. 


15 


Parabolic  extrapolation  is  done  in  the  following  way.  Suppose  the  imaginary  parts  of 
three  “complexified”  k-values  are  and  k2  and  the  corresponding  RCS  values  for  one 

observation  angle  are  ao,  cti,  and  <72,  respectively.  We  compute  the  Lagrange  interpolation 
polynomial  through  the  points  (fco,^o),  (&i,cti),  and  (k2,cr2)  and  set  k  equal  to  zero  in  this 
quadratic  polynomial  in  k( i.e.,  we  extrapolate  to  the  real  axis)  obtaining 

ff(0)  =  po<?o  +  P\a\  +  P202  (29) 


where  p0  =  kik2/[(k0  -  fci)(&0  -  ^2)],  Pi  =  k0k2/[(ki  -  fco)(*’i  ~  *2)],  and  p2  =  k0ki/[(k2  - 
k0)(k2  ~  /ci )].,  Parabolic  extrapolation  is  second-order  accurate  provided  the  function  cr(k) 
is  “sufficiently  smooth”  (has  three  continuous  derivatives).  This  means  that  for  uniform 
spacing  (A k  =  k2  —  fcj  =  fcj  —  k0 )  the  error  in  a  due  to  the  extrapolation  is  0((Afc)2). 
Linear  extrapolation  is  only  0(Ak )  if  a  has  two  continuous  derivatives  in  k.  In  practice, 
we  have  found  that  the  accuracy  of  the  extrapolation  procedure  is  even  better  than  these 
estimates  would  suggest.: 

In  Fig.  2,  we  show  some  RCS  plots  for  scattering  from  a  two-dimensional  conducting 
circular  cylinder  for  linear  extrapolation  from  “complexified”  k  values.  To  remove  the 
effect  of  the  “complexification”  in  computing  the  RCS  in  the  far  field,  wc  must  rescale 
FTVff)  where  rpp  is  the  value  of  r  in  the  far  field,  typically  about  10,000.  If  k  = 
krt  +  ik,„,  where  kun  is  the  imaginary  term  added  to  k  to  accelerate  convergence  of  the 
iterative  procedure,  then  we  rescale  E3zcat(rpp )  by  multiplying  it  by  exp( ktm rpp)  from  the 
asymptotic  expansion  of  H^\ki„,rpp).  The  RCS  is  then  computed  using  equation  (27) 
and  the  rescaled  values  of  E3cat(rpp).  Extrapolation  of  the  RCS  to  the  real  axis  is  also  done 
using  the  rescaled  values.  In  Fig.  2,  the  linearly  extrapolated  values  (from  k  =  80  +  0.5 i 
and  k  —  SO  +  1.0?)  of  the  RCS  give  exact  comparisons  with  the  unextrapolated  ones  (for 
k  =  80).  The  parabolically  extrapolated  values  give  the  same  excellent  agreement.  The 
number  of  iterations  required  to  converge  to  an  error  tolerance  of  10-3  for  these  three  cases 
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are  79,  26,  and  16,  respectively,  and  again  for  faster  convergence  it  pays  to  “complexify”. 
Note  also  that  large  condition  numbers  not  only  imply  slow  convergence  but  also  poor 
accuracy,  giving  another  reason  to  “complexify”.  Of  course,  one  could  always  solve  the 
CFIE  rather  than  equation  (4)  and  avoid  the  problem  of  large  condition  numbers  due 
to  resonance.  Unfortunately,  solving  the  CFIE  for  TE  polarization  requires  dealing  with 
the  second  derivative  of  the  free  space  Green’s  fuction,  which  has  a  nasty  singularity. 
“Complexification”  is  easier! 

In  Fig.  3a  and  Fig.  3b,  we  plot  the  RCS  for  scattering  from  a  conducting  elliptical 
cylinder  using  linear  extrapolation  (a  =  2, b  =  1,  k  =  100  +  0.5?,  k  =  100  +  1.0?)  with  an 
incident  plane  wave  at  90°.  As  can  be  seen  in  Table  1,  n  =  1600  for  these  cases.  The 
agreement  between  extrapolated  values  and  unextrapolated  values  of  the  RCS  is  exact. 
Again  we  merely  show  the  extrapolated  values  in  the  figures  to  avoid  too  much  congestion. 
Parabolic  extrapolation  using  k  =  100  +  0.3?,  100  +  0.6?,  100  +  0.9?  gives  the  same  excellent 
results. 

As  we  have  stated,  the  FMM  algorithm  of  this  paper  is  0(??4/3 )  per  conjugate  gradient 
iteration.  This  means  that  the  CPU  time  is  proportional  to  ??4/3,  i.e., 

T  =  avAX  n4^3  (30) 

where  we  have  emphasized  that  the  proportionality  constant  a  is  a  fuction  of  the  computer 
being  used  (in  this  case  a  VAX  6410  in  double  precision).  Using  equation  (30)  and  some 
of  the  data,  in  Table  1,  we  can  obtain  a  rough  estimate  for  the  value  of  a \/.\x-  In  this  case 

qvax  ~  -003  (31) 

Suppose  now  we  would  like  to  estimate  the  CPU  time  required  on  our  VAX  to  com¬ 
pute  the  electromagnetic  scattering  problem  for  an  example  having  10.000  unknowns.  We 
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assume  that  without  too  much  “complexification”  convergence  occurs  in  10  iterations  or 
less.  Thus,  for  this  example 

Ttot  «  ovAx(10)n4/3  =  .03(10, 000)4/3  =  1.8  hours  (32) 

Although  this  last  CPU  time  may  seem  large,  we  are  only  using  a  VAX  in  double  precision, 
and  a  vectorized  version  of  this  code  on  a  CRAY  computer  would  do  the  same  calculation  in 
minutes.  Furthermore,  ordinary  Gaussian  elimination  for  a  matrix  of  order  10,000  would 
require  about  (10, 000)3/3  operations  and  would  use  ponsiderably  more  CPU  time  than 
that  in  equation  (32). 

Finally,  the  most  important  statement  about  the  FMM  and  our  approach  to  “com¬ 
plexification”  and  extrapolation  of  the  RCS  values  for  complex  ks  to  the  real  axis  is  that 
for  complex  k  ( Re(k )  >  0  and  Im(k)  >  0),  the  solution  to  equation  (4),  a  second  kind 
integral  equation,  is  unique  for  sufficiently  smooth  scatterers  and  consequently,  the  con¬ 
dition  number  of  A  is  bounded14  and  independent  of  n.  This  means,  when  we  use  our 
extrapolation  procedure  (either  linear  or  parabolic),  Nlter  is  small  (say,  Nticr  <  10)  and 
independent  of  n  for  error  tolerances  of  10-3.  Of  course,  if  the  accuracy  requirements  are 
increased,  so  will  Ar,<er,  but  in  either  case  the  FMM  is  globally  0(n4/3)  and  we  can  write 

T  °  —  Niter  Ot computer  n4^3  (33) 

where  Niter  is  only  a  function  of  e,  the  error  tolerance,  and  not  a  function  of  n  ,  and 
a computer  is  a  function  of  the  computer  speed  for  arithmetic  operations  and,  also  inde¬ 
pendent  of  n  .  For  further  proof,  we  present  more  evidence  in  Table  2,  where  we  have 
set  the  “complexification”  to  5 i  in  all  cases  and  maintained  approximately  10  points  per 
wavelength  for  various  geometries.  Note  the  small  number  of  iterations  to  convergence, 
indicating  a  bounded  condition  number  independent  of  n.  Finally,  in  Table  3,  we  list  some 
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preliminary  results  of  the  FMM  code  on  the  CRAY  2.  This  code  is  not  yet  fully  vectorized, 
so  we  don’t  give  CPU  times. 

Although  we  have  only  considered  circular  and  elliptical  cylinders  for  the  test  cases 
in  this  paper,  the  code  can  handle  any  closed  two-dimensional  metal  scatterer  having  a 
unique  outward  normal  at  every  point  on  C.  Therefore,  at  present,  we  cannot  deal  with 
the  trailing  edge  of  an  airfoil  that  comes  to  a  point  (non-existence  of  dG/dn  at  the  trailing 
edge);  currently,  we  must  round  this  region  off.  However,  we  are  working  on  replacing  the 
current  at  this  singular  point  with  its  asymptotic  expansion,  allowing  us  then  to  deal  with 
such  points.  We  are  also  considering  thin  coatings  about  metal  scattcrers  using  higher 
order  impedance  boundary  conditions,  similar  to  our  treatment  in  Ref. [15]. 

Our  closing  observation  is  that  in  computational  electromagnetic  scattering,  almost 
anyone  can  obtain  reasonable  results  using  a  bad  algorithm  if  ka  is  small.  It  is  much 
more  difficult  to  get  good  results  for  large  values  of  ka,  and  it  is  doubly  difficult  to  do 
so  efficiently-  We  hope  that  by  choosing  moderate  values  of  ka,  we  have  demonstrated 
the  robustness  of  the  FMM  algorithm  for  effectively  solving  electromagnetic  scattering 
problems. 
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Table  1 

Fast  Multipole  Results  for  Scattering  from  Two-Dimensional  Conducting 
Circular  and  Elliptical  Cylinders 
(DEC  VAX  6410  Computer,  Double  Precision) 


kb 

a 

b 

n 

Niter 

CPU 

Err 

50 

1 

1 

500 

10 

66 

11.99 

.441  (-3) 

80 

1 

1 

800 

10 

79 

22.35 

.856  (-3) 

80+  1.6i 

1 

1 

800 

10 

13 

22.43 

.200  (-3) 

80  +  2.4i 

1 

1 

800 

10 

9 

21.73 

.655  (-3) 

80  +  3.2i 

1 

1 

800 

10 

8 

21.89 

.198  (-3) 

80  +  4.8i 

1 

1 

800 

10 

6 

21.27 

.750  (-3) 

80  +  6.4i 

1 

1 

800 

10 

6 

21.22 

.223  (-3) 

150  +  6.4i 

1 

1 

1500 

10 

6 

51.94 

.633  (-3) 

100 

2 

1 

1600 

10.12 

184 

34.89 

.992  (-3) 

100 +  2i 

2 

1 

1600 

10.12 

11 

35.57 

.293  (-3) 

100 +  3i 

2 

1 

1600 

10.12 

9 

36.53 

.478  (-3) 

100 +  4i 

2 

1 

1600 

10.12 

8 

35.84 

.494  (-3) 

100  +  6i 

2 

1 

1600 

10.12 

6 

35.35 

.998  (-3) 

50 

3 

1 

1200 

10.73 

95 

32.61 

.950  (-3) 

50  +  2i 

3 

1 

1200 

10.73 

10 

33.98 

.863  (-3) 

100 +  2i 

3 

1 

2237 

10 

12 

86.93 

.639  (-3) 

k  =  Wavenumber 
a  =  Semimajor  axis 
b  =  Semiminor  axis 

n  =  Total  number  of  unknowns  (sample  points  on  scatterer  boundary) 
nx,  =  Number  of  unknowns  per  wavelength  of  incident  radiation 
Njter  =  Number  of  iterations  to  convergence 
CPU  =  CPU  time  per  iteration,  in  seconds 
Err  =  Relative  Error 


NOTATION:  .441  (-3)  means  0.441  X  10'3 


Table  2 

Fast  Multipole  Method:  CPU  Time  for  Fixed  “Complexification” 

( Scattering  from  Two-Dimensional  Conducting  Circular  and  Elliptical 

Cylinders) 

(DEC  VAX  6410  Computer,  Double  Precision) 


kb 

a 

b 

n 

Niter 

CPU 

Err 

50  +  5i 

1 

1 

500 

6 

11.97 

.251  (-3) 

100 +  5i 

1 

1 

1000 

6 

22.46 

.973  (-3) 

150  +  5i 

1 

1 

1500 

7 

42.28 

.313  (-3) 

200  +  5i 

1 

1 

2000 

7 

51.90 

.686  (-3) 

50  +  5i 

2 

1 

791 

7 

22.22 

.998  (-3) 

150 +  5i 

2 

1 

2372 

8 

78.59 

.408  (-3) 

100 +  5i 

3 

1 

2237 

9 

86.28 

.999  (-3) 

k  =  Wavenumber 
a  =  Semimajor  axis 
b  =  Semiminor  axis 

n  =  Total  number  of  unknowns  (sample  points  on  scatterer  boundary) 
Niter =  Number  of  iterations  to  convergence 
CPU  =  CPU  time  per  iteration,  in  seconds 
Err  =  Relative  Error 

NOTATION:  .251  (-3)  means  0.251  X  10’3 


Table  3 

Fast  Multipole  Method:  CRAY  Results  with  Fixed  “Comnlexification” 
(Scattering  from  Two-Dimensional  Conducting  Circular  a.:u  Elliptical 

Cylinders) 

(CRAY-2  Computer,  Single  Precision) 


kb 

a 

b 

n 

Niter 

Err 

200 +  5i 

1 

1 

2000 

7 

.686  (-3) 

300 +  5i 

1 

1 

3000 

8 

.331  (-3) 

400 +  5i 

1 

1 

4000 

8 

.803  (-3) 

500 +  5i 

1 

1 

5000 

9 

.327  (-3) 

100 +  5i 

2 

1 

1582 

7 

.997  (-3) 

200  +  5i 

2 

1 

3163 

6 

.995  (-3) 

300 +  5i 

z. 

1 

4744 

6 

.991  (-3) 

100 +  5i 

3 

1 

2237 

8 

.995  (-3) 

200  +  5i 

3 

1 

4473 

7 

.997  (-3) 

100 +  5i 

4 

1 

2916 

9 

.998  (-3) 

50  +  5i 

10 

1 

3554 

8 

.996  (-3) 

k  =  Wavenumber 
a  =  Semimajor  axis 
b  -  Semiminor  axis 

n  =  Total  number  of  unknowns  (sample  points  on  scatterer  boundary) 
Niter  ~  Number  of  iterations  to  convergence 
Err  =  Relative  Error 


NOTATION:  .686  (-3)  means  0.686  X  10*3 


FIGURE  CAPTIONS 


Fig.  la.  Radar  cross  section  (RCS)  for  a  circular  cylinder  of  unit  radius.  Incident  field  is 
a  plane  wave  with  wave  number  k  =  80.  Solid  line  shows  solution  from  the  FMM  code, 
while  dotted  line  is  that  from  the  series  solution. 

Fig.  lb.  RCS  of  previous  case  from  the  FMM  code  for  more  observation  angles. 

Fig.  2.  RCS  for  a  circular  cylinder  of  unit  radius.  Incident  fields  are  plane  waves  with 
with  wave  numbers  k  =  80,  k  =  80  +  .5 i,  and  k  =  80  +  l.Oi.  Solid  line  shows  solution  from 
the  FMM  code  with  no  extrapolation  ( k  —  80).  Dotted  line  shows  solution  of  the  linearly 
extrapolated  RCS  from  k  =  80  +  .5 i  and  k  =  80  +  l.Oi. 

Fig.  3a.  RCS  results  from  the  FMM  code  for  an  elliptical  cylinder  with  semimajor  axis 
a  —  2  and  semiminor  axis  6=1  linearly  extrapolated  from  incident  field  plane  waves  at 
90°  with  wave  numbers  k  =  100  4-  .5i  and  k  =  100  +  l.Oi. 

Fig.  3b.  RCS  of  previous  elliptical  case  for  more  observation  angles. 
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ABSTRACT 


Higher  order  impedance  boundary  conditions  for  thin  coatings  about  closed  conduc¬ 
tors  in  two  dimensions  are  derived  using  Fourier  integral  techniques.  Using  a  single-layer 
potential  and  these  impedance  boundary  conditions,  second-kind  weakly  singular  integral 
equations  are  derived  that  model  TE  electromagnetic  scattering  problems.  These  integral 
equations  are  solved  using  Nystrom’s  method  and  approximately  fourth-order  convergent 
quadrature  formulas. 


♦ 
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I.  INTRODUCTION 

Consider  a  two-dimensional  closed  perfect  electrical  conductor  coated  with  a  thin  layer 
of  dielectric  and/or  magnetic  material.  The  classical  way  of  solving  the  electromagnetic 
scattering  problem  from  such  an  object  is  to  develop  an  integral  equation  in  which  the 
contour  of  integration  contains  both  the  conductor  and  the  outer  surface  of  the  dielectric. 
The  difficulty  with  this  approach  is  that  as  the  thickness  of  the  dielectric  layer  approaches 
zero,  an  ill-conditioned  equation  may  result.  In  addition,  the  size  of  the  discrete  linear 
system  is  twice  as  large  as  in  the  method  described  below.  Our  procedure  will  translate 
the  boundary  condition  on  the  surface  of  the  conductor  to  the  dielectric-air  interface  by 
developing  an  impedance  boundary  condition  on  the  interface.  The  resulting  integral 
equation  will  only  have  to  be  integrated  along  the  interface,  thereby  reducing  the  number 
of  unknowns  for  the  discrete  problem  by  a  factor  of  two  and  possibly  removing  the  ill- 
conditioning  caused  by  grid  points  on  the  conductor  and  dielectric  being  too  close  together. 
Other  work  in  this  area  can  be  found  in  Rojas  and  Al-hekail1,  Senior  and  Volakis2’3, 
Harrington  and  Mautz4,  Karp  and  Kara!5,  Senior6,7,8,  and  Barkeshli  and  Volakis9. 
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II.  DERIVATION  OF  HIGHER  ORDER  IMPEDANCE  BOUNDARY  CONDITIONS 


Consider  Helmholtz’s  equation  written  in  co-ordinates  normal  (n)  and  tangential  (t) 
to  the  scatterer,  i.e., 


d2utot  .  d2utot  .  tot 


dn2  +  dt 2 


+  k{utot  =  0 


(1) 


where  fci  is  the  wave  number  in  the  dielectric  (&i  =  kairy/e Here  the  superscript  tot 
denotes  the  total  field.  Thus 

utot  =u  +  uinc  (2) 

where  inc  denotes  the  incident  field.  No  superscript  is  used  for  the  scattered  field.  The 
boundary  condition  on  the  conductor  for  the  TE  polarization  case  is 


dutot 

dn 


(3) 


and  on  the  dielectric-air  interface  we  have 


du?0*  _  £2  du\ot 
dn  e\  dn 


(4) 


where  ei  is  the  dielectric  constant  in  the  thin  layer  and  62  is  that  in  air.  We  assume  the 
layer  thickness  is  8. 

Taking  Fourier  integral  transforms  around  the  scatterer  and  assuming  periodic  condi¬ 
tions,  we  can  derive  the  following  equation,  where  “  denotes  Fourier  transformation: 


d2utot 
dn 2 


=  (a;2  -  k2)utot  s  autot,  0  <n<8 


(5) 


Let  Uq0(  denote  the  unknown  value  of  utot  at  n  =  0  (on  the  metal-dielectric  interface). 
Then  we  obtain 

utot(  0)  =  ulot  (6a) 


4 


(66) 


* 


k 


a2rulo<(o) 

dn*r 


=  aruiot 


dutoi(  0)  ...  ^^(O)  __  A  r_19 
an  an2--+i  ’  ’ 


(6c) 


where  we  have  used  equations  (3)  and  (5).  Expanding  in  a  Taylor  series  in  6,  we  have 


utot(8  -  0)  =  u*6t  + 


82au  o°‘ 


(7) 


auto<(6  -  o) 

dn 


=  8auf,ot  + 


83a2U(,ot 

— r^  + 


Substituting  (8)  into  (4)  gives  (at  n  =  8  +  0) 


(8) 


du^/dn  =  e28(aut0ot  +  (a£)2Uo0</6  +  •  •  ‘)/ei 

=  e28aut20t/e1  +  0(82)  (9) 

where  we  have  used  equation  (7)  and  the  continuity  of  utot.  Note  that  the  subscript  2 
represents  the  point  n  =  6  +  0.  Taking  inverse  transforms  of  equation  (9)  yields  the  desired 
impedance  boundary  condition 

=  -6'2(d*ur/dt>  +  kluy) 

dn  d  '  '  ^ 


In  terms  of  scattered  and  incident  fields  equation  (10)  may  be  re-written  as 


du2  ,  ci  \^2u2  .  .9  •<  duinc 


dn 


+ «(•/*)  Fs£ + -  *(«./«.)  + *W] 


at2 


=  9 


(ii) 


where  we  have  dropped  the  0(82)  term  in  equation  (11). 
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Now  introduce  the  single-layer  potential 


«(*)  =  J  $(*>  y)<f>(y)  ds(y)  (12) 

where  x  and  y  are  two-dimensional  points,  C  represents  the  outer  contour  (around  the 
dielectric),  and  $  is  the  two-dimensional  free  stream  Green’s  function  given  by 

$(z,  y )  =  iH^l)(k \x  -  y|)/4  (13) 

Here  denotes  the  first  kind  Hankel  function  of  order  zero,  and  k  =  kair.  If  equa¬ 
tion  (12)  is  substituted  into  (11)  and  the  appropriate  jump  conditions10  are  enforced  at 
the  dielectric  surface,  we  can  derive  the  following  weakly  singular  second-kind  integral 
equation: 

<f>(x)-2  f  [ d$(x,y)/dn(x)  +  6(e2/€ii)(d2$(x,y)/dt2(x)  +  kl$(x,y)j\<f>(y)ds(y ) 

J  c 

=  -2 g(x)  (14) 

See  Colton  and  Kress10  for  more  details  on  derivations  of  integral  equations  in  the  form 
of  equation  (14).  If  S  is  set  to  zero,  equation  (14)  reduces  to  the  TE  polarization  case 
for  scattering  from  a  perfect  electrical  conductor.  If  the  term  d2$/dt2  is  removed  from 
equation  (14),  we  have  the  standard  first  order  impedance  boundary  condition  integral 
equation10.  Finally,  as  written,  we  call  equation  (14)  the  second  order  (0(S2))  impedance 
boundary  condition  integral  equation  for  modeling  a  thin  coating  about  a  metal  conductor. 
The  second  derivative  term  obviously  models  the  curvature  of  the  scatterer. 

The  above  derivation  can  easily  be  extended  to  multiple  dielectric  surfaces.  Suppose 
we  have  m  overlapping  coatings  with  physical  parameters  (Si,  e;,  k{)  (i  =  1,2, .'. .  ,m).  Let 
6  =  Si  +  62  -I - b  Sm  and  set  =  e,+i/ej  where  em+i  =  eair-  Then  following  the  same 
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steps  as  above,  we  can  derive  the  following  impedance  boundary  condition  at  the  m-th 
dielectric-air  interface: 

m 

du/dn  =  i(d2u/dt2  +  k2ju)  +  g  +  0(62)  (15) 

i= i 

where  we  have  dropped  the  subscript  2,  and  g  is  now  defined  to  be  equal  to  the  sum  in 
equation  (15)  for  the  incident  field.  An  integral  equation  analogous  to  equation  (14)  can 
be  easily  written. 

In  addition,  fourth  and  higher  order  impedance  boundary  conditions  may  be  developed 
by  allowing  more  terms  in  the  Taylor  series  expansions  (7)  and  (8).  For  example,  the  fourth 
order  impedance  boundary  condition  is  given  by 

du/dn  +  ( e2/e1)(63/6)[d2/dt 2  +  k\\[d2u/dt2  +  Jfejtt]  =  g  (16) 

where  —  g  now  takes  the  form  of  the  second  term  in  equation  (16)  for  the  incident  field  and 
the  subscript  2  has  been  dropped. 
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III.  NUMERICAL  METHOD 


We  discretize  equation  (14)  using  Nystr  om’s  method11  and  the  approximately  fourth- 
order  convergent  quadrature  formulas  that  handle  logarithmic  singularities  derived  in 
Refs.[12,13,14].  The  main  advantages  of  Nystrom’s  method  over  the  method  of  moments15 
or  finite  elements  are  that  matrix  fill  is  less  expensive  and  higher  order  convergent  quadra¬ 
ture  formulas  are  easier  to  employ  than  corresponding  high  order  basis  functions  such  as 
piecewise  cubic  Hermite  polynomials  or  cubic  B-splines16. 

The  resulting  linear  system  is  then  solved  by  Gaussian  elimination  with  partial  piv¬ 
oting  or  generalized  conjugate  residual(GOR)17  iterative  methods.  The  latter  algorithm 
performs  well  for  moderate  condition  numbers  (/c  <  1000).  For  larger  condition  numbers, 
we  use  the  theory  of  ’’complexification”.  See  below  and  Ref.  [14]  where  we  have  seen  that 
this  theory  is  applicable  even  at  resonance  frequencies  (k  — »  oo). 
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IV.  COMPUTATIONAL  RESULTS 


Assume  the  incident  field  is  in  the  form  of  a  plane  wave  with  the  incident  angle  0,  i.e., 

u,nc(x )  =  ex\>[ik(xi  cos 0  +  xi  sin/3)]  (17) 

where  x  =  (x\,xi).  For  our  first  example  we  consider  a  thin  dielectric  around  a  metal 
circular  cylinder,  since  the  analytic  solution18  is  available  for  this  case. 
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The  “resonance  problem"  is  that  at  certain  values  of  the  wave  number  k  (the  resonant  k's),  the 
second-kind  integral  equation  for  solving  scattering  problems  can  become  extremely  ill- 
conditioned.  This  adversely  affects  both  the  accuracy  and  speed  of  numerical  solutions.  We 
consider  transverse-magnetic  scattering  from  a  conductor  (Dirichlet  problem).  The  integral 
equation  (derived  using  double-layer  potentials)  is  discretized  using  approximately  fourth- 
order  convergent  quadrature  formulas.  At  resonant  k ’s  for  circular  and  elliptical  scatterers,  we 
find  very  large  condition  numbers  for  the  discrete  matrices  [up  to  0(  107)  ],  generally  leading 
to  poor  solutions.  We  apply  two  approaches  to  alleviate  the  resonance  problem.  The  first  is  to 
use  a  different  integral  equation,  based  on  both  single-  and  double-layer  potentials.  This  leads 
to  low  condition  numbers  and  good  solutions  at  resonant  k.  The  second  method  is  to  use  the 
original  second-kind  integral  equation,  introduce  a  small  imaginary  part  in  k,  and  extrapolate 
back  to  the  real  axis.  Solutions  obtained  by  the  two  methods  are  in  excellent  agreement.  The 
extrapolation  technique  will  be  particularly  useful  in  the  case  of  the  exterior  Neumann 
problem,  when  the  application  of  the  first  technique  will  be  numerically  more  difficult.  By 
solving  the  resonance  problem,  we  ensure  that  fast  and  accurate  solutions  are  obtainable  at  any 
arbitrary  wave  number. 


I.  INTRODUCTION 

A.  What  is  the  resonance  problem? 

Briefly,  the  resonance  problem  is  that  at  certain  values 
of  wave  number,  the  second-kind  integral  equation  for  solv¬ 
ing  scattering  problems  can  become  extremely  ill-condition¬ 
ed.  The  presence  of  resonances  adversely  affects  both  the 
accuracy  and  the  speed  of  the  solution  method:  The  accura¬ 
cy  is  affected  because  of  the  ill-conditioning  of  the  problem. 
The  speed  may  be  affected  because  large  condition  numbers 
tend  to  make  it  difficult  for  iterative  solution  techniques  to 
converge. 1  The  resonance  problem  must  be  solved  if  reliable 
and  efficient  solutions  are  to  be  obtainable  at  any  arbitrary 
frequency. 

The  resonance  problem  arises  in  many  numerical  inte¬ 
gral-equation  methods  for  solving  scattering  problems.  Pre¬ 
vious  investigators  have  discussed  the  problem  for  moment- 
method  solutions  of  the  magnetic-field  integral  equation 
(MFIE),2  electric-field  integral  equation  (EFIE),'  and 
combined-field  integral  equation  (CFIE).24 

In  this  paper  we  concentrate  on  alleviating  the  reso¬ 
nance  problem  in  the  Nystrom-method  solution  of  the  sec¬ 
ond-kind  integral  equation  for  electromagnetic  scattering 
derived  using  double-layer  potentials.  The  basic  numerical 
method  is  described  by  Murphy,  Rokhlin,  and  Vassiliou5 
and  is  summarized  in  Secs.  I  B  and  II  A  below. 

B.  A  more  precise  definition 

Consider  a  two-dimensional  closed  conducting  object  in 
the  {x ,y)  plane.  We  treat  transverse-magnetic  (TM )  scatter¬ 
ing  from  such  an  object.  By  TM  we  mean  that  the  electric 


field  of  the  incident  radiation  is  in  the  z  direction:  i.e.,  paral¬ 
lel  to  the  axis  of  the  scatterer  (in  transverse  electric,  or  TE, 
scattering  the  magnetic  field  is  parallel  to  the  axis). 

One  way  of  modeling  TM  scattering  from  a  closed  con¬ 
ducting  object  defined  by  the  curve  C  is  to  introduce  the 
double-layer  potential 

«(-*)  =  f  *(y)ds(y).  (1) 

Jc  dv(y) 

where  x  and  y  are  points  in  the  plane,  x  is  outside  C  and  y  is 
on  C.  <t>  is  the  two-dimensional  Green's  function 

<!>(■*, y)  H I,"  (k  \x  —  y; )  (2) 

4 

where  k  is  the  wave  number  of  the  incident  radiation  and 
H'vu  is  the  Hankcl  function  of  the  first  kind  of  order  zero. 
Physically,  u  represents  the  scattered  electric  field  E.  per¬ 
pendicular  to  the  plane  containing  the  scatterer.  u  satisfies 
Helmholtz’s  equation,  with  a  radiation  condition  at  infinity 
and  with  the  boundary  condition  determined  by  the  incident 
electric  field  £ 

u=/U)  =  onC.  (3) 

This  is  a  Dirichlet  problem.  Colton  and  Kress'1  show  that  the 
function  'P  satisfies  the  integral  equation 

V(a)  +  2  f  V(y)ds(y)  =  2 f(x),  *eC,  (4) 

Jc  dv{y) 

Unfortunately,  Eq.  (4)  does  not  have  a  unique  solution  for 
all  values  of  k.  In  fact,  nonunique  solutions  exist  for  those 
values  of  k  where  the  interior  Neumann  problem 
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VJu  +  k  7p  =  0  in  C, 

(5a) 

i2.-0  one 

(5b) 

dv 

has  a  nontrivial  solution.  These  values  of  k  are  referred  to  as 
interior  resonances  of  the  Dirichlet  (TM)  problem,  or  inte¬ 
rior  Dirichlet  eigenvalues. 

Numerically,  when  one  tries  to  solve  Eq.  (4)  at  or  near  a 
resonance,  the  condition  number  k  of  the  resulting  discre¬ 
tized  linear  system  becomes  large.  The  solution  becomes 
much  more  sensitive  to  computer  round-off  errors,  and  iter¬ 
ative  solution  procedures  have  more  difficulty  converging. 
In  short,  efficient  solution  techniques  and  accurate  results 
become  much  more  difficult  to  achieve. 

C.  What  can  be  done? 

Below,  we  study  the  condition  number  of  the  discrete 
system  for  the  case  when  C  is  a  circle  and  the  case  when  C  is 
an  ellipse.  We  then  outline  two  approaches  to  alleviate  the 
resonance  problem:  ( 1 )  using  a  combined-potential  equa¬ 
tion  in  place  of  Co  (4)  and  (2)  introducing  a  small  imagi¬ 
nary  part  in  k  and  extrapolating  back  to  the  real  axis. 

II.  THE  ILL-CONDITIONING  OF  SCATTERING 
PROBLEMS  NEAR  A  RESONANCE 

A.  Numerical  method:  Discretization  and  quadrature 
formulas 

We  discretize  Eq.  (4)  using  Nystrom’s  method;  that  is, 
we  approximate  the  integral  directly  with  a  quadrature  for¬ 
mula.  We  use  roughly  fourth-order-convergent  quadrature 
formulas  tl.at  handle  logarithmic  singularities.5,7  The  high- 
order  accuracy  is  important  because  in  the  vicinity  of  a  reso¬ 
nance,  solutions  are  extremely  sensitive  to  round-off  and 
truncation  errors. 

Roughly  fourth-order-convergent  quadrature  formulas 
that  handle  logarithmic  singularities  at  x  —  0  are  developed7 
by  starting  with  the  Euler-Maclaurin  formula  with  the  sin¬ 
gular  point  removed.  That  is,  let 

I fix)dx  =  *  (,?/<'■  >  -^T1)  -  TT*  ™ 

(6) 

where  h  —  l/n  and  x(  =  i/n.  To  correct  for  the  singular 
point,  a  concentration  of  points  of  the  form  2,6=  ,  AJ(.Xj )  is 
introduced  in  the  first  interval,  where  Xj  —  ft  A>  for 
j  =  1,2,...  ,6.  The  derivative  term  is  approximated  by  the 
one-sided  difference  formula 

/'(*„)=—[  3/(x„ )  -  4/(x„  _,)+/(*„_  2 1 J  -  ( 7 ) 

In 

Combining  these  terms  yields 
f  f{x)dx  =  h  (  £  f{x, )  +  £  A./(Xj )  -  ) 

JO  \/  ~  I  /  a  I  l 

-  (/(*„ -  2 )  -  4/(x„ .  , )  +  3 f(Xll ) }). 

(8) 

The  unknown  weights  Aj  (J  =  1,2,...  ,6)  are  determined  by 


solving  the  6x6  linear  system  that  results  when  Eq.  (8)  is 
assumed  exact  for  the  following  candidate  functions  /(x): 
l,x,x\logx,xlogx,  and  xMogx,  using  analytic  integration 
rules.  Once  computed,  the  quadrature  weights  can  be  stored 
and  looked  up  numerically  when  needed.  The  proof  of  ap¬ 
proximately  fourth-order  convergence  can  be  found  in 
Rokhlin.7 

8.  The  condition  number 

The  discrete  system  is  of  the  form 

Ax  =  b,  (9) 

where  A  is  a  nonsymmetric  complex  matrix  of  order  n  and  x 
and  b  are  vectors  of  length  n.  The  condition  number  k  is 
defined  as 

*  =  II  IM  "‘II.  OO) 

where  ||  |j  can  be  any  finite-dimensional  norm.  We  choose 
the  2-norm,  in  which  case  we  have8 

«•  = - ,  (11) 

^min 

where  o-mM  and  amm  are,  respectively,  the  largest  and  small¬ 
est  singular  values  of  the  matrix^.  We  can  thus  compute  the 
condition  number  of  the  matrix  by  performing  a  singular 
value  decomposition.  This  is  a  somewhat  expensive  oper¬ 
ation  to  perform  routinely,  but  it  is  indispensible  for  the  pres¬ 
ent  study.  As  a  resonance  frequency  is  approached,  crmm  ap¬ 
proaches  zero,  and  the  condition  number  approaches 
infinity. 

C.  Resonances  of  a  circular  conductor  of  unit  radius 

For  a  circle  of  unit  radius,  the  resonant  wave  numbers 
are  determined  from  the  roots  of  the  derivatives  of  Bessel 
functions;  that  is, 

J'n(km)=  0,  n  =  0,1,2 .  m  =  0,1,2 .  (12) 

In  Table  I,  we  show  some  computational  results  for  k  with  k 
at  or  near  a  resonance,  with  the  scatterer  being  a  circle  of  unit 
radius.  The  value  of  n  denotes  the  number  of  sample  points 
along  the  boundary,  separated  by  equal  arc  lengths.  The 
number  of  points  per  wavelength  is  n/k.  The  three  resonant 
wave  numbers  k ,,  k2,  and  k}  in  Table  I  have  values 
3.831  705  970  2,  7.015  586  669  8,  and  19.615  858  510  5,  re¬ 
spectively.  More  roots  of  J can  be  found  in  Abramowitz 
and  Stegun.1’ 

Notice  in  Table  I  that  a  slight  movement  of  k  into  the 
complex  plane  reduces  the  condition  number  k  by  orders  of 
magnitude,  even  for  the  highest  wave  number  kt.  We  will 
use  this  fact  in  our  extrapolation  technique. 

D.  Resonances  of  an  elliptical  conductor 

For  an  ellipse,  the  resonant  wave  numbers  are  deter¬ 
mined  from  the  roots  of  the  derivative  of  a  modified  Mathieu 
function.  We  generated  these  roots  ourselves:  to  our  knowl¬ 
edge  they  are  not  satisfactorily  tabulated.  Mathieu  functions 
arise  often  in  problems  with  elliptical  geometr>  and  have  an 
extensively  developed  theory,'0  which  we  cannot  present 
here.  Since,  however,  they  are  not  as  generally  known  or  as 
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TABLE !.  Condition  numbers*  for*  circle  of  radius  a  =  1,  using  Eq.  (4). 
Notes:  0.10(3)means0.1xl0\  n  is  the  number  ofsample  points  on  the 
scatterer  boundary,  k  is  the  wave  number,  and  a  is  the  radius;  a  =  I. 


Number  of  points 
Wavelength 


n 

ka 

K 

on  scatterer  boundary 

40 

k, 

0.34790(5) 

10.4 

80 

k, 

0.92251  (6) 

20.8 

too 

A, 

0.29271  (7) 

26.1 

200 

k, 

0.50210(7) 

52.2 

100 

k,  +  0.05  / 

0.20355  (2) 

26.1 

100 

k,  +  0.025/ 

0.19506  (2) 

26.1 

too 

k,  +  0.01/ 

0.10163  (3) 

26.1 

too 

k,  +  0.002/ 

0.50803  (3) 

26.1 

100 

k,  +  0.001/ 

0.10160(4) 

26.1 

too 

A,  +  0.0004/ 

0.25401  (4) 

26.1 

too 

k,  +  0.0002/ 

0.50804  (4) 

26.1 

40 

0.38677  (4) 

5.7 

80 

ki 

0.96858  (5) 

lt.4 

100 

ki 

0.27361  (6) 

13.7 

200 

K 2 

0.83871  (7) 

27.4 

too 

kj  +  0.05/ 

0.20126(2) 

13  7 

100 

k}  +0.001/ 

0.10062  (4) 

13.7 

100 

A, 

0.63729  (4) 

5.1 

200 

*, 

0.16200  (  6) 

10.2 

200 

k,  +  0.05/ 

0.20194(2) 

10.2 

widely  tabulated  as  Bessel  functions,  we  discuss  them  brief¬ 
ly.  When  we  apply  the  method  of  separation  of  variables  to 
Helmholtz’s  equation  in  elliptical  coordinates,  we  obtain  the 
following  equations: 

d*v 

— ~+  (A  —  2qcos2u)y=z0  (13a) 

du 2 

and 

^-rr  —  (A  —  2?cosh  2u)y  =  0,  (13b) 

du1 

which  are,  respectively,  the  Mathieu  and  modified  Mathieu 
equations.  In  the  present  problem,  the  equation  of  interest  is 
the  second  one,  the  modified  Mathieu  equation.  Solutions  to 
this  equation  are  the  modified  Mathieu  functions,  which 
have  the  form" 


Ce2„  (u,q)  =  V  A  jJn,cosh  2 ru, 

(14a) 

r  *  O 

ee 

Ce2n  +  ,  ( u,q )  =  £  A  ^"++,  "cosh(2r+  l)u, 

(14b) 

r*  O 

oc 

Se2„+ ,  (a,<?)  -  X  ^2J"+V,s>nh(2r+ l)w, 

r  «  0 

(14c) 

Se2n  +  2(«.?)  =5  X  £2?++22,sinh(2r-f  2)u. 

( 14d) 

r  _  0 


The  notation  Ce  and  Se  (on  which  there  are  variations) 
arose  from  the  notion  of  “elliptic  cosine"  and  “elliptic  sine.” 

The  coefficients  A  and  B  must  be  found  through  recur¬ 
rence  relationships.  First,  however,  one  must  determine  the 
“characteristic  number”  A  for  which  a  periodic  solution  ex¬ 
ists.  One  must  do  this  for  each  value  of  the  parameter  q  (giv¬ 
en  u,  which  is  fixed  by  the  eccentricity  of  the  ellipse  in  ques¬ 


tion).  To  find  the  characteristic  number  one  must  find  the 
root  of  a  continued-fraction  equation.  Our  approach  is  to 
obtain  tabulated  values  of  the  characteristic  numbers,10  in¬ 
terpolate  them  using  a  natural  cubic  spline,  and  use  the  inter¬ 
polated  value  as  an  initial  guess  in  the  method  of  bisection  to 
find  the  root  of  the  continued-fraction  equation.  This  gives 
us  the  characteristic  number  for  any  value  of  q  to  high  preci¬ 
sion.  We  can  then  determine  the  coefficients  and  calculate 
the  values  of  the  derivative  of  the  modified  Mathieu  function 
for  various  values  of  q,  and  use  the  secant  method  to  obtain 
the  derivative’s  root.  The  root  q0  is  related  to  the  resonant 
wave  number  k  via  q0  =  (c2/4)  A: J,  where  2 c  is  the  interfocal 
distance  for  the  ellipse  in  question.  We  have  generated  roots 
for  the  function  Ce £  for  use  in  this  study. 

Our  present  computer  implementation  requires  that  the 
sample  points  on  the  boundary  of  the  scatterer  be  equispaced 
in  arc  length.  In  the  case  of  a  circular  scatterer,  this  is  simple: 
As  =  aA0,  where  As  is  the  arc  length  interval  and  a  is  the 
radius  of  the  circle.  However,  for  other  closed  contours  a 
numerical  algorithm  must  be  invoked  to  prxluce  equally 
spaced  points.  This  resampling  procedure,  which  uses  fast 
Fourier  transforms,  perturbs  the  discrete  problem  very 
slightly  (typically  0.01%  or  less);  however,  these  small  per¬ 
turbations  are  often  great  enough  so  that  the  resulting  linear 
system  has  a  different  resonance  from  the  original  problem. 

To  obtain  the  resonant  wave  numbers  of  the  discrete 
problem  resulting  from  a  slightly  perturbed  scatterer  shape, 
we  use  the  secant  method  to  drive  the  minimum  singular 
value  <rmn  to  0.  Suppose  the  resonant  wave  number  is  de¬ 
fined  as  the  root  of  a  nonlinear  function  <7min  =f(k).  Then, 
using  the  secant  method,  an  improved  root  is  given  by 


£  t»  +  i )  __  £  m>  _  jj-i") 


k -  k 


(n-  I) 


T("> 


(15) 


where  the  superscript  (ft)  denotes  the  iteration  counter.  In 
practice,  we  have  found  that  only  one  or  two  iterations  are 
required  to  obtain  an  accurate  resonant  wave  number  for  an 
ellipse,  when  we  use  the  root  of  the  derivative  of  *he  appro¬ 
priate  modified  Mathieu  function  to  obtain  a  starting  value 
for  the  algorithm.  Table  II  lists  the  continuous  and  discrete 
resonant  wave  numbers  for  a  number  of  ellipses.  The  values 
k4,  A's,  'k6,  and  k7  refer  to  the  discrete  k ’s  in  this  table  in 
order  from  top  to  bottom,  and  are  used  in  the  tables  follow¬ 
ing. 

Table  III  lists  some  condition  numbers  for  various  val¬ 
ues  of  k  near  a  discrete  resonance  for  an  ellipse  with  semima¬ 
jor  axis  a  and  semiminor  axis  b.  We  compute  the  number  of 
points  per  wavelength  using  an  approximate  formula  for  the 
perimeter  of  an  ellipse,  which  yields 


TABLE  H.  Resonant  wave  numbers  for  various  ellipses.  Notes:  a  is  the 
semimajor  axis,  b  the  semiminor  axis,  and  k  the  wave  number. 


Continuous  kb 

Discrete  kb 

a 

b 

3.380  123  647  98 

3.377  841  278  90  (=  A.) 

2 

1 

3.467  850  747  85 

3.467  387  362  51  (=  A,) 

1.5 

1 

3.562  999  874  37 

3.558  639  339  65  (  =  k„ ) 

1.25 

1 

3.301  504  871  94 

3.300  352  24101  (  =A,) 

3 

1 
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TABLE  111.  Condition  numbers  k  for  various  ellipses,  using  Eq  (4). 
Notes.  0.10  (3)  means  0.1  X  10'.  //  is  the  number  of  sample  points  on  the 
scatterer  boundary,  A  is  the  wave  number,  a  is  the  semimajor  axis,  and  b  is 
the  scniiminor  axis. 


n 

kb 

a 

b 

K 

Number  of  points 
Wavelength 
on  scatterer  boundary 

150 

*4 

2 

1 

0.15107  (8) 

28.07 

150 

k4  -4*  0,002/ 

2 

1 

0.50537  (3) 

28.07 

150 

k,  +  0.001/ 

2 

1 

0.10107  (4) 

28.07 

150 

k,  +  0.02/ 

2 

1 

0.50542  (2) 

28.07 

150 

A,  +0.01/ 

2 

1 

0.10108  (3) 

28.07 

150 

k,  +  0.2/ 

2 

1 

0.51228  (1) 

28.07 

150 

A\,  +  0.1/ 

2 

1 

0.10141  (2) 

28.07 

150 

a. 

1.5 

1 

0.12572  (8) 

33.93 

150 

A,  +0.002/ 

1.5 

1 

0.50698  (3) 

33.93 

150 

A,  +  0.001/ 

1.5 

1 

0.10139  (4) 

33.93 

150 

A* 

1.25 

1 

0.36088  (7) 

37.24 

150 

A„  +  0.002, 

1.25 

1 

0.50202  (3) 

37.24 

150 

A„  +0.001/ 

1.25 

1 

0.10040(4) 

37.24 

150 

^7 

3 

1 

0.19492(8) 

20.32 

150 

A,  +0.002/ 

3 

1 

0.52592  (3) 

20.32 

150 

A,  +0.001/ 

3 

1 

0.10522  (4) 

20.32 

points  _ n _ 

X  kjiia'  +  b2)  ' 

Table  III  shows,  as  does  Table  I,  the  bad  conditioning 
that  occurs  near  resonance.  Again,  notice  how  a  slight  move¬ 
ment  of  k  into  the  complex  plane  reduces  k  by  orders  of  mag¬ 
nitude. 

III.  SOLVING  THE  RESONANCE  PROBLEM  USING  THE 
COMBINED-POTENTIAL  EQUATION 

A  different  integral  equation  from  Eq.  (4)  can  be 
formed  by  combining  single-  and  double-layer  potentials. 
Let 

«(*)  =  f  (  -  iy<t>(x,y)  W)cfrO>),  (17) 

Jc  \  dv(y)  / 

where  775*0  is  an  arbitrary  real  number  such  that  . 


77  Re  k> 0. 

(18) 

We  use  77  =  l  in  our  computations.  Using  (17),  Colton  and 
Kress6  derive  the  integral  equation 

V  +  KV  -  irjSW  =  2 f 

(19) 

where  K  and  S  are  operators  given  by 

(*¥)(*)  =  2  |  'i'(y)ds(y), 

Jc  dv(y) 

(20) 

(SYMx)  =  2  £<I>(Jc,y)'I'(y)dy(y). 

(21) 

Equation  ( 19)  is  the  modified  Dirichlet  problem  (TM)  that 
has  solutions  for  all  wave  numbers  k  satisfying  Im  k> 0* 
Equation  (19)  is  analogous  to  the  combined-field  integral 
equation  (CFIE)  of  electrical  engineering  practice,  but  it  is 
not  the  same.  We  refer  to  it  in  this  paper  as  the  combined- 
potential  equation. 


TAHLE  IV.  Condition  numbers*  for  11  circle  of  radius  u  =  I  using  ihv  com¬ 
bined-potential  equation  |Eq.  (19)  |.  Notes:  //  is  the  number  of  sample 
points  on  the  scatterer  boundary,  A  is  die  wave  number,  and  a  is  the  radius; 
o=l. 


11 

ka 

K 

Number  of  points 
Wavelength 
on  scatterer  boundary 

40 

A, 

3.8933 

10.4 

80 

A, 

3.8930 

20.8 

100 

A, 

3.8930 

26.1 

40 

A, 

7.0648 

5.7 

80 

A, 

7.0507 

11.4 

100 

*1 

7.0505 

13.7 

200 

A, 

7.0505 

27.4 

If  we  discretize  Eq.  ( 19 )  using  the  same  highly  accurate 
quadrature  formulas  we  used  for  Eq.  (4),  a  linear  system 
similar  to  (9)  results,  but  with  a  much  better-conditioned 
coefficient  matrix  A.  We  illustrate  this  point  in  Tables  IV 
and  V.  Equation  (19),  the  combined-potential  equation,  is 
somewhat  more  complicated  than  Eq.  (4 ) ,  the  integral  equa¬ 
tion  obtained  by  using  only  double-layer  potentials.  How¬ 
ever,  the  linear  system  is  so  much  better  conditioned  that  the 
additional  work  required  to  discretize  the  integral  equation 
is  justified.  Since  most  iterative  methods  for  solving  the  lin¬ 
ear  system  (9 )  have  a  convergence  time  that  is  an  increasing 
function  of/c,  the  condition  number,  k  should  be  minimized. 

Unfortunately,  in  the  transverse-electric  (TE)  case, 
which  is  a  Neumann  problem,  the  combined-potential  equa¬ 
tion  analogous  to  Eq.  ( 19)  contains  the  second  derivative  of 
the  Green’s  function  <t>.  This  makes  the  equation  much  more 
difficult  to  solve  numerically.  Thus,  we  investigate  an  alter¬ 
native  approach  to  solving  the  resonance  problem. 

IV.  SOLVING  THE  RESONANCE  PROBLEM  BY 
EXTRAPOLATION  FROM  THE  COMPLEX  PLANE 

In  this  simple  but  powerful  approach,  we  take  the  reso¬ 
nant  wave  number,  add  to  it  a  small  imaginary  part,  and 
solve  the  scattering  problem.  We  repeat  this  for  a  different 
small  imaginary  part.  We  then  have  two  different  solutions 
that  we  can  use  to  extrapolate  to  a  solution  corresponding  to 
the  resonant  wave  number  with  zero  imaginary  part.  We 
show  an  example  in  Fig.  I.  Here,  we  use  solutions  for 
k  —  -+  /0.002  and  k  =  kt  +  z'0.001  to  obtain  the  radar 

TABLE  V.  Condition  numbers  k  for  various  ellipses  using  the  combined- 
potential  equation  [Esq  (19)].  Notes’  n  is  the  number  of  sample  points  on 
the  scatterer  boundary,  A  is  the  wavenumber,  01s  the  semimajor ar is,  and  b 
is  the  semiminor  axis. 


n 

kb 

a 

b 

K 

Number  of  points 
Wavelength 
on  scatterer  boundary 

150 

A, 

2 

1 

3.2821 

28.07 

150 

A, 

1.5 

1 

3.3914 

33.93 

150 

As 

1.25 

1 

3.4160 

37.24 

150 

A, 

3 

1 

3.5451 

20.32 
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converge  rapidly  for  such  small  values  of  /c.121' 

The  RCS  results  for  the  ellipse  a  =  2,  b  =  1  at  resonance 
k,  are  plotted  in  Fig.  2  [this  curve  is  the  same  on  a  plot  of  this 
size  whether  computed  via  Eq.  (19)  or  via  the  extrapolation 
technique).  These  results  are  compared  in  the  figure  with  the 
RCS  results  computed  at  resonance  using  Eq.  (4)  with  no 
extrapolation.  The  enormous  difference  between  the  two 
curves  highlights  the  possibly  disastrous  results  that  can  be 
obtained  when  condition  numbers  reach  values  of  0(10  ). 
Blindly  solving  Eq.  (4)  without  considering  resonance  is  not 
wise. 


FIG.  1.  Radar  cross  section  (RCS)  for  a  circle  of  unit  radius  at  resonant 
wave  number  k, .  Solid  line  shows  solution  obtained  using  the  combined- 
potential  equation  ( 19).  Dotted  line  shows  solution  obtained  using  extrapo¬ 
lation  from  the  complex  plane.  The  error  between  the  two  methods  is  less 
than  0.002  dB  and  cannot  be  seen  on  the  plots.  The  series  solution  for  a  circle 
is  also  indistinguishable  on  the  plot. 

cross  sc-tion  (RCS)  for  a  circle  of  unit  radius  at  resonant 
wave  number  kx .  As  can  be  seen  in  Fig.  1,  the  agreement 
between  this  method  and  that  using  Eq.  (19)  described 
above  is  excellent.  The  error  between  the  two  methods  is  less 
than  0.002  dB  and  cannot  be  seen  on  the  plots. 

As  another  example,  we  consider  an  elliptical  scatterer 
with  parameters  a  *=  2,  b  =  1  ( where  a  is  the  semimajor  axis 
and  b  is  the  semiminor  axis),  n  =  150,  and  k  =  kA,  illumi¬ 
nated  by  a  plane  wave  with  incidence  angle  90°  (broadside 
incidence).  Three  separate  linear  extrapolations  are  com¬ 
puted:  ( 1 )  k4  +  /0.002  and  kA  +  /0.001,  (2)  &4  +  /0.02  and 
kK  -f  iO.Ol,  and  (3)  kA  +  /0.2and  kA  +  /0. 1 .  The  maximum 
RCS  errors  ((or  100  observation  angles)  between  the  solu¬ 
tions  obtained  using  these  three  extrapolations  and  those  ob¬ 
tained  using  the  combined-potential  equation  (19)  are,  re¬ 
spectively,  0.003,  0.002,  and  0.146  dB.  Even  the  last 
extrapolation,  a  crude  one  that  probably  would  not  be  used 
in  practice,  yields  a  good  solution  at  resonance.  Note  also 
from  Table  III  that  the  condition  numbers  k  for  the  cases 
kt  +  i0.2  and  k<  -f  iO.  1  are  almost  as  small  as  that  for  kA  in 
Table  V.  Similarly  good  results  are  obtained  for  the  other 
ellipses  (with  different  values  of  a  and  b)  in  Table  III.  Itera¬ 
tive  methods  such  as  the  generalized  conjugate  residual 
(GCR)  method  for  solving  linear  systems  can  be  expected  to 


FIG.  2.  RCS  results  for  #n  ellipse  with  semimajor  axisa  =  2and  semiminor 
axis  6  =  1  at  resonance  k, .  The  extrapolation  technique  (dotted  line)  and 
combined-potential  equation  technique  agree  very  well  with  each  other  and 
would  be  indistinguishable  if  different  symbols  were  not  used.  The  other 
solid  line  represents  the  incorrect  results  computed  at  resonance  using  Eq 
(4)  with  no  extrapolation. 


V.  SUMMARY 

Wc  have  solved  the  problem  of  resonance  in  integral- 
equation  scattering  methods  using  two  different  approaches: 
(1)  By  solving  a  combined-potential  equation  [Eq.  (19)) 
and  (2)  by  extrapolation  from  the  complex  plane.  In  this 
paper  we  have  solved  a  Dirichlet  (TM)  problem  for  which 
the  combined-potential  equation  is  relatively  easy  to  treat 
numerically,  so  that  the  extrapolation  technique  may  not 
appear  to  offer  any  overwhelming  advantages.  However, 
when  we  solve  the  Neumann  (TE)  problem,  the  combined- 
potential  equation  analogous  to  Eq.  ( 19)  contains  the  sec¬ 
ond  derivative  of  the  Green’s  function  d>.  This  makes  the 
equation  much  more  difficult  to  solve  numerically.  In  the  TE 
case,  the  extrapolation  technique  will  offer  a  significant  ad¬ 
vantage. 

By  solving  this  problem,  we  have  assured  the  possibility 
of  obtaining  accurate  and  efficient  solutions  at  all  frequen¬ 
cies.  We  have  thus  significantly  strengthened  a  highly  accu¬ 
rate  implementation  of  the  second-kind  integral  equation  for 
scattering. 
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NUMERICAL 

SECOND-KIND-INTEGRAL-EQUATION 
SOLUTIONS  OF  ELECTROMAGNETIC 
SCATTERING  PROBLEMS 


Indexing  terms:  Electromagnetic  tvatici,  Scattering,  Numerical 
methods  and  number  theory.  Radar  cross-sections 

A  new,  highly  accurate  numerical  method  based  on  second- 
kind  integral  equations  has  been  developed  to  solve  electro¬ 
magnetic  scattering  problems  for  closed  conducting  bodies  in 
two  dimensions.  The  method  is  approximately  fourth-order 
convergent,  owing  to  the  use  of  accurate  new  quadrature 
formulas 


Introduction :  We  have  developed  a  second-kind  integral  equa¬ 
tion  solver  (the  SKIE  method)  for  transverse-magnetic  (TM) 
electromagnetic  scattering  from  perfect  electrical  conductors 
of  arbitrary  closed  geometry  in  two  dimensions.  The  method 
uses  accurate,  roughly  fourth-ordcr-convergcnt  quadrature 
formulas.  The  resulting  discrete  matrix  has  a  condition 
number  bounded  by  a  constant  as  the  sampling  is  refined  (for 
a  given  scatterer  size  and  nonresonant  frequency). 

Method:  We  begin  with  a  second-kind  integral  equation 
derived  from  Helmholtz’s  equation 

{  eG^dn  ~  ^  dq  +  0) 


We  solve  a  Dirichlet  problem,  with  <f>  the  prescribed  value  of 
the  potential  on  the  closed  boundary  B  of  the  scatterer.  In  the 
TM  case  we  are  solving,  4>  is  the  electric  field  normal  to  the 
plane  of  the  scatterer,  q  and  p  are,  respectively,  source  and 
field  evaluation  points  on  B,  dG/dn  is  the  normal  derivative  of 
the  Green’s  function  for  Helmholtz’s  equation  and  p  is  the 
unknown  double-layer  distribution  on  the  boundary  of  the 
scatterer.  Once  p  is  known,  the  field  at  a  point  exterior  to  the 
scatterer  can  be  obtained  by  evaluating 


<t>(p)  = 


'  dOjt,  q) 

R  dtl 


p(q)  dq 


(2) 


Eqn.  1  is  discretiscd  by  replacing  the  integral  with  quadrature 
formulas,  leading  to  a  matrix  equation  that  can  be  solved 
numerically.  Wc  have  solved  the  linear  system  using  both 
Gaussian  elimination  and  an  iterative  generalised  conjugate 
residual  technique.1  Roughly  fourth-order-convergent  quadra¬ 
ture  formulas  that  handle  logarithmic  singularities  at  x  =  0 
are  develop'  /  by  starting  with  the  Euler-Maclaurin  formula 
with  the  si1  ,u!ar  point  removed.  That  is,  let 

i 

J/W  dx  =  *[  £  /(*,)  -  ~]  -  ^  *7  '(X.)  (3) 

0 

where  h  =  1/n  and  x,  =  i/n.  To  correct  for  the  singular  point,  a 
concentration  of  points  of  the  form  Yj- 1  f(Xj)  is  introduced 
in  the  first  interval,  where  Xj  =  j7i/6  for  j  =  1,  2,  ....  6.  The 
derivative  term  is  approximated  by  the  one-sided  difference 
formula 


/'<*„)  =  Yh  [3/W  _  «> 


(4) 


Combining  these  terms  yields 

i 


I  m  dx  =  h\i  f(x,)  +  i  ).jf(Xj)  -fix.)/ 2 
J  (i-i  y  i 


i  [/(*„- 2)~  4/(x,.,)  +  3/(x,)]).  (5) 


\ 


The  unknown  weights  Ajj  =  1,  2,  6)  are  determined  by 

solving  the  6  x  6  linear  system  that  results  when  eqn.  5  is 
presumed  exact  for  the  following  candidate  functions  /(x):  1, 
x,  x2,  log  :t,  x  log  x,  and  x2  log  x,  using  analytic  integration 
rules.  Once  computed,  the  quadrature  weights  can  be  stored 
and  looked  up  numerically  when  needed.  The  proof  of 
approximately  fourth-order  convergence  can  be  found  in 
Reference  2. 

Our  method  involves  a  direct  discretisation  of  the  integral 
using  a  quadrature  formula,  and  does  not  explicitly  involve  an 
expansion  in  basis  functions  as  in  the  method  of  moments.5 
Common  implementations  of  the  method  of  moments  typi¬ 
cally  use  pulse,  triangle  or  low-order  trigonometric  basis  func¬ 
tions.  Such  choices  of  basis  functions  lead  at  best  to  first  or 
second-order  accuracy,  provided  the  singularity  in  the  Green’s 
function  is  handled  adquately.  (In  theory,  it  is  possible  to 
obtain  higher  accuracy  by  using  higher-order  basis  functions.) 

Results:  The  method  as  implemented  at  present  can  calculate 
TM  electromagnetic  scattering  from  closed,  perfectly  conduct¬ 
ing  bodies  in  two  dimension.  Some  bistatic  RCS  Circulations 
are  shown  in  Figs.  1  and  2.  Fig.  1  shows  the  results  for  a 
cylinder  of  size  ka  =  20  where  a  is  the  radius,  and  Fig.  2 
shows  the  results  for  an  ellipse  of  size  ka  -  27-7  where  a  is  the 
semimajor  axis. 

Accuracy:  One  way  to  verify  the  accuracy  of  the  numerical 
method  is  to  compare  numerical  results  with  known  analytical 
solutions  for  canonical  shapes,  as  done  for  the  cylinder  of 
Fig.  1.  We  also  adopt  an  approach  where  the  accuracy  can  be 
verified  for  any  shape  of  object.  The  approach  ;s  essentially 
one  of  testing  how  closely  the  method  verifies  Green’s  second 
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Fig.  I  Bislanc  RCS  ( TM )  for  conducting  cylinder  of  size  ka  =  20, 
where  a  is  radius 

To  show  detail,  lobe  al  0°  (forward  scattering)  is  cu.  off  in  plot,  and 
plot  shows  only  first  half  of  azimuth  range 
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identity.  Suppose  that  we  arc  solving  an  exterior  Dirichlet 
problem  on  some  two-dimensional  boundary,  and  the  right- 
hand  side  of  the  problem  is  equal  to  the  field  generated  by  a 
dipole  source  located  inside  the  boundary.  The  solution  of  this 
Dirichlet  problem  at  receivers  exterior  to  the  boundary  must 
be  equal  to  the  field  of  the  source  at  those  recr  -s,  as  if  the 
boundary  were  not  present.  A  relative  error  can  oc  computed 
between  the  fields  calculated  by  solving  the  Dirichlet  problem 
and  the  fields  of  the  source. 


Fig.  I  Bistatic  RCS  {TM)  for  broadside  incidence  on  conducting  ellipse 
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FI*  3  Relative  accuracy  against  fineness  of  sampling  for  ellipse  (no* 
same  ellipse  as  in  Fig.  2) 

The  accuracy  of  the  method  has  been  verified  for  two- 
dimensional  scatterers  of  many  different  shapes.  An  example 
of  the  results  of  one  accuracy  test,  for  an  ellipse  of  eccentricity 
0-87  and  ka  =  5  (a  is  the  semimajor  axis)  is  shown  in  Fig.  3. 
The  average  value,  over  360°,  of  the  absolute  value  of  the 
relative  error  as  defined  above  is  plotted  against  the  number 
of  points  per  wavelength  on  the  scatterer.  As  the  number  of 


Fig.  4  Condition  number  of  discrete  matrix  against  number  of  sample 
points  per  wavelength  ellipse  of  Fig  3 


points  per  wavelength  is  increased  from  5  to  60,  the  relative 
error  decreases  from  order  10* 1  to  order  I0"‘,  The  dotted 
line  shows  fourth-order  convergence,  and  is  drawn  for  refer¬ 
ence.  At  10  points  per  wavelength,  the  accuracy  is  of  order 
t0'3. 

Condition  number:  Fig.  4  illustrates  another  important  pro¬ 
perly  of  the  SKIE  method.  For  the  ellipse  mentioned  above  in 
connection  with  Fig.  3,  Fig.  4  plots  the  condition  number  of 
the  matrix  against  the  number  of  sample  points  per  wave¬ 
length.  As  Fig.  4  shows,  the  condition  number  is  low,  and  it 
does  not  change  as  the  number  of  sample  points  is  increased. 
The  SKIE  method  produces  problems  that  are  well- 
conditioned  away  from  resonances.  A  low  condition  number 
is  advantageous  both  because  it  indicates  a  stable  method4 
and  because  it  tends  to  lead  to  faster  convergence  for  iterative 
methods  for  linear  systems.3  Note  that  the  condition  number 
is  bounded  by  a  constant  for  a  given  scatterer  size  and  non- 
resonant  frequency. 
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